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TECHNICAL MEMORANDUM 

A NUMERICAL METHOD FOR INTERFACE PROBLEMS 
IN ELASTODYNAMICS 

I. INTRODUCTION 


In ground motion studies of earthquakes, the principal interest is in the effect 
on a localized region of the Earth. This region can take the form of a structure 
foundation, a region behaving nonlinearly near a structure, or a region that has 
significantly different soil properties from its surroundings. The difficulty in dealing 
with this type of problem analytically, is the infinite nature of the Earth, in relation 
to the region of interest, and its energy absorbing characteristics. In order to 
address these problems, either very large finite element models are required or 
absorbing boundaries must be used. Large finite element models are computationally 
inefficient and standard absorbing boundaries often do not accurately represent the 
energy dissipation of the Earth. 

This problem can be defined in terms of an interior finite region and an 
exterior infinite region. This type of problem is called an "Interface Problem." A 
method, being investigated recently, to study these problems combines finite element 
and boundary integral methods [1,3]. The exterior region is handled by means of a 
boundary integral formulation and the interior is handled with finite element methods. 
Continuity conditions are placed across the boundary involving the displacement field 
and the corresponding tractions. 

In this study, this technique is applied to a specific problem. The goal is to 
develop a computer implementation of the method and determine its accuracy and con- 
vergence. The problem chosen is that of an anti-plane SH shear wave incident upon 
a semi-cylindrical alluvial valley. A closed form analytical solution for this problem 
is available [2] and will be the basis for the comparisons made. 

The specific problem presented here is a simple one. The study, however, has 
two significant goals: (1) to compare the approximate method with exact analytical 

results, and (2) to establish a firm baseline for further extensions of the method. 


II. MATHEMATICAL FORMULATION 


A. Problem Formulation 

In order to formulate the problem in question, the region of interest must first 
be defined. Consider the geometry of Figure 1» The region ft + is to represent 
homogeneous elastic material with mass density p + and Lame’s constants X + , p + . The 
region ft represents a possibly inhomogeneous cylinder with properties p - , X~ , and p 
depending on the spatial coordinates x^ and Xg. C is a traction-free surface and r 
is the interface between ft and ft + . 



A displacement field must be defined next. The ease of antiplane strain due to 
an incoming antiplane shear wave shall be studied. While this represents a rather 
restricted type of excitation, the corresponding response will have qualitative similari- 
ties with that due to more general types of plane wave excitation. This approach can 
then be useful in gaining physical insight into the problem of evaluating the effects 
of local geology on ground motion. 

The single non-vanishing, antiplane displacement is assumed to be time-periodic 
of the form : 

W(x r x 2 ,t) s; Re[w t (x 1 ,x 2 ) e lwt ] (1) 


Also, w°(x 1 ,x 2 ) corresponds to the driving field. w° will consist of an incoming wave 

defined in the half space x 2 > 0 and its reflection from the surface x 2 - 0, taken to 

be traction free. For specifying w° it will be assumed that the entire halfspace is 
homogeneous, with properties equal to those of ft + , That is, it will be assumed that 
the excitation displacement field, or free-field displacement, w° is known in the 
absence of the obstacle ft. 

For convenience, a displacement w is defined such that: 
t 

! w in ft 

t ( 2 > 

w -w° in ft + 


Thus, in ft, w represents the total displacement; that is, the displacement due to the 
free-field excitation, w°, and the effects of the obstacle ft. In ft + , w represents only 
the effects of the obstacle ft without the free-field excitation. In other words, in ft + , 
w represents the scattered field only. 


The problem then is to find the time periodic displacement field w, inside and 
outside the obstacle ft, such that it satisfies the three following* conditions. First, 
since both regions ft and ft + are linearly elastic, the displacement field w must satisfy 
a generalized, reduced wave equation in both of these regions. Second, w must 
satisfy several boundary conditions; namely, the displacements and tractions must be 
continuous across r and C is traction free. And third, w must satisfy a radiation 
condition in ft + . These conditions can be expressed as follows: 


(y~ W ) + (y~ w ) + p 

X 1 X 1 x 2 x 2 

“ w 2 w = 0 

in 

ft 

(3) 

y + V 2 W + p + 

w = 0 


in 

+ 

ft 

(4) 

w t+ - w 1 " = 0 ; 

+ t + 
y w n - 

y” wl" - 0 
n 

on 

r 

(5,6) 

y W n = ° 



on 

c 

(7) 

w ^ r" 1/2 e" ir 

as r = x. 2 

+ x 2 *+• « 



(8) 


where 


a> = frequency of free-field motion, w° 

t+ 4- 

w = total displacement on ft T side of boundary r 
w 1 '" = total displacement on ft side of boundary r 

w^ f - normal derivative of total displacement on ft + side of boundary F 
w^~ = normal derivative of total displacement on ft side of boundary r . 


Taking equation (2) into consideration, equations (5) and (6) can be rewritten as: 


w“ = w + w° 


— - + + + 0 
V w n = y w n + v w n 


(9) 

( 10 ) 


respectively, where: 

w n ° = normal derivative of free-field displacement. 
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Equations (3), (4), (7), (8), (10), and (11) define a well-posed problem for w in Q 
and £2 + . 


It is desirable, from a computational point of view, to reformulate the problem 
in a way that eliminates the need to deal with an infinite exterior region, which is 
both time and space consuming. This is done at the expense of introducing two 
boundary functions on r. It begins with some potential theory for the Helmholtz 
equation? 

V 2 v + p + w 2 v = 0 , X 2 > 0 (11) 

subject to the radiation condition (8) and the traction-free boundary condition, 


v x = 0 on x 2 ~ 0 


G is defined by the formula: 


G(x,y) = r | H o ( 2) (klrsD + h 0 (2) (12) 

where 

k 2 = to 2 p + /jj + 

y* = (y r -y 2 ) 

H q 2(.) = the Hankel function of the second kind and zeroth order. 

Physically, G can be interpreted as the displacement produced at x due to a steady- 
state point load of frequency wand unit amplitude at y. 

It is assumed that the solution v in can be represented in terms of a simple 
layer with density $ ; that is , it is assumed that a function 0 exists , defined over r , 
such that : 

v(x) = f G(x,y) <J>(y) dSy in . (13) 

r 


Thus, the displacement v at every point within can be represented in terms of a 
density function <p defined only on r. It can be shown [ 5] that G is a continuous 
function and, therefore, the limiting displacement v* on r is given by 

v± (x) = f G(x,y) <p(£) dSy on T . (14) 

r 
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In addition, under the assumption that r is smooth, the limiting normal derivative is 
[5] 


+ i C 8G(x,jr) 

v n “(x) s ± - 2 - $(x) + J $(jr) dSy on r 


(15) 


r x 


The plus and minus denote limits from S2 + and Si, respectively, and n is the unit 

■X 

normal vector on r at point x (Fig. 1). Then equations (4), (5), and (6) can be 
replaced in the original problem statement by a combination of equations (9), (10), 
(14), and (15). This yields 



G(x,y) <Ky) dSy + w° 


on r 



1 + 

■2 » 


4>(x) + y + 


/ 


9G(x,y) 



<J>(y) dSy + + w n ° 


on r 


In order to arrive later at a symmetric discretized formulation of the problem, 
the following analogues of the Helmholtz formulas are introduced. If v satisfies equa- 
tions (11) and (12) in Q then, 


f 3G(x,y) f 

1/2 v” = J — g n - y ~ v~(y) dSy - J G(x,u) v n "(y) dSy in SI . (16) 


If v satisfies equations (11) and (8) in Sl + then, 


+ f + f 3G(x,y) 

1/2 v = J G(x,y) v n (y) dSy - J v (y) dSy in fl 


(17) 


By equations (9) and (10) w + = w"~w° and w + = y“/y + w "-w°, respectively. By 
substituting these results into equation (17), n 


1/2 w“ + 


/ 


9G(x,y) 

9ny 


w 


(y) dSy - / G(x,y)^w n "^ (y) dSy 


f f 3G(x,y) 

" 1/2 w ° - J G(x,y) w n °(y) dSy + J w°(y) dS Z » 0-8) 


is obtained. Now w° is a solution of equation (11) in S2, hence by equation (16), 
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(• 3G(x,jr) f 

1/2 w ° = J --5— — w°( j;) dSjr - J G(x,jr_ w n °(}r) dS X . 


So the right hand side of equation (18) becomes w° and can be rewritten, 


f 3 G(x > 5 r) 

1/2 W ~ + J -3HT - w " 


(X) <IS£ - J G(X >X ) ^ w n ") (jr) dSjr 


W° . (19) 


This leads to the following problem. Find (w, y” w n ~, <j>) such that 


(y~ w ) + <y" w ) + p~ w 2 w * 0 in ft 

x 3, x i x 2 2 


= / G (x, 


y) <Ky) + w ° on r 


(20a) 


y" w n~ = 1 ^ 2 ^ K x ) 


x) + y + / 


9G(x,£) 


<Kx) dSy f u + w n ° on r (20c) 


y“ w = 0 on C 
n 


1/2 w ~ + f 

r 


9G(x,y) 


■<Z> dSZ - / GCx.jr) U n "j ( X ) ds z = 


(20e) 


The only set of points now used are those in ft and r . The set of points in ft is 
completely described by points in r by the use of the Green’s function. It has been 
shown in Reference 1 that this problem has, in general, a unique solution provided 
certain critical values of w are excluded . 

B . Variational Principle 

Now standard Galerkins' ideas are used to obtain a variational formulation 
which will be suitable for discretization by the finite element method. First, multiply 
equation (20a) by a test function <5w and integrate over ft using the divergence 
theorem and equation (20d). The result is, 


J {y“( w 6w + w 6 w ) - w p~ w 5 w } dx + I 
ft X 1 x l x 2 x 2 r 


+ / »" v 


<5w dS = 0 


(21a) 
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Next, multiply equation (20b) by a test function <5 y“ w n ~ and integrate over r; 


([{ 


G(x,y) <Ky) dSy - w"(x) + w°(x) 6 y“ w n "(x) dx 




(21b) 


Perform similar calculations for equations (20c) and '(20c) and combine the resulting 
equations linearly with (21a) and (21b) to obtain: 


- f {y“(w v 6 w + w v <5 w v ) - w 2 p“ w <5 w) dx + / y“ w “ <$w dS 

L J -‘ 1 1 2 2 " r J 

- 1/2 uu G(x,y) <p(x) dSy - w~(x) + w°(x)| <5 y~ w n ~(x) dS j 
+ 1/2 £/* |l/2 y + <j>(x) + y + / 

+ 1/2 |l/2 y* w"(x) + y + f — g--- w"(y) dSy - f G(x,y) (y“ w n “)(y) dSy 

” y ' w° | 6 tf>(x) dS^ a o 


0G(x,y) + 

3n' x ~" ' W dS X + U w n ° - y V*l 5 w " ds 


The problem is now in a form suitable for finite element approximation. This 
is done by defining the various functions in equation (22) using finite element shape 
functions as follows: 


I w 

w(x) = [N(x)] T {w} = [N^(x) T : N r (z) T ] < J 1 

\ r / 

Hx) = [Q(X)J T {0 


y" w n “(x) = A(x) = IQ(x)] {A} . 


(23a) 

(23b) 

(23c) 


The corresponding test functions are defined using the same shape functions. 
Using these shape functions equation (22) leads to the matrix equation: 
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K 


RO 


K 


nr 


K 


rn 


K 


0 

0 


rr 

A 


A J 

0 


B 


B 


C 


T 



W R 


o , 


V 


r r 

, 

j 

A 1 


' h 


* 

. 

V’ 


( 24 ) 


where 


[K] = - f {[N 3 v“[N ) T + [N ] y~[N ] T - w 2 [N] T p“[N)} dx 

J n X 1 1 2 x 2 

[A3 - 1/2 / [Q] [N r ) T dS 


[B3 “ 1/2 y 
[CJ - - 1/2 


1/2 / [Q] [N r ] A dS + 

r r 'r r 


/* f (" 3G(x,y)"l m 

Jj h^J l V ds * ds 


iT 


/ / [Q3 [G(X,X>1 CQ1* dSjr dS 

r r 


{f } - - 1/2 y + / [N r 3 w n ° dS 

* i* 


{f A > = 1/2 jT CQ] w° dS 


{f^> = 1/2 y 


/ 


[Q] w° dS 


The complete development of equation (24) from equation (22) is shown in Appendix 

A. 


For convenience U> and {$} are now condensed out, so the equation is in terms 
of {w} only. 


K 


RR 


K 


Rr 


K r r K r r + D ] 


M 


1 i 



W R 

ii 

0 


w r 


F 

r 

, 


(25) 


is obtained where 


[D rr ] =-[A) T [CJ 1 [B] - [BJ T [C T ] _1 [A] 


{F r > = (f r > - [A] T IC]” 1 {fo - [B] T [C T r 1 {f A > . 


S 

t 

J 


.ri 


, V 


From equation (24), it can be seen that [K] is essentially the stiffness matrix 
of the cylinder; that is, the force at some point in the cylinder due to a ur u t dis- 
placement at some other point fe, the cylinder, [CJ represents the discretized com- 
pliance of the halfspace. [A] and [B] are matrices which couple the cylinder 
boundary and the halfspace boundary and guarantee the continuity of displacements 
and tractions across r . 

From equation (25), [D rr l is the force at a point on the boundary due to a 

force or traction applied at some other point on the boundary allowing for the* 
presence of the halfspace medium. r r } is essentially a representation t>£ a nonlocal 

absorbing boundary developed from a Green’s function rather than a spring and dash- 
pot model. 


III. APPLICATIONS 


The primary objective of this study was to develop a computer code to imple- 
ment the previously described formulation and to check its accuracy. In order to do 
this, a baseline reference was needed, The most commonly cited problem of this type 
deals with a semicircular deposit with homogeneous material properties, for which an 
exact solution is available [21. The program developed, therefore, used these 
characteristics. Extensions to more general problems can be made once the method 
has been checked. 

The free-field antiplane shear wave, w°, was taken as 
-i(wt - kXj^ sin e Q ) 

tv 0 (x^,x 2 ) - e cos(kx 2 cos 0 Q ) (26) 

where 


0 O - angle of incidence 

, 2 _ 2 + , + 
k - a) p /y 


Figure 2 illustrates x^, x 2 , and 0 q and their sign conventions. 
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With the geometry being circular, polar coordinates were used for their con- 
venience. By substituting polar coordinates and simplifying, the incident wave 
becomes; 


, . i kr sin(0+(L) -i kr sin (0-0^) 

w° = e" 1 “ t 1/2 [o 0 + o 0 ] 


(27) 


and the normal derivative; 


8 w° _ 9 W° _ iwt 

9 n “ 3 r 


i kr sin (0+0 ) 

1/2 (i k sin(0+0 o ) e 


i kr sin(o-0 ) 
- i k sin( 0- 8 q ) e ° 3 


The Green’s function (12) can be expressed as 


(28) 


G(x,y) = \ H q (2) [2 kr [sin |] 


where 


(29) 


AO = angular distance between x and y. 

4 . 

Note that this form of the Green's function satisfies the wave equation in n and the 
radiation condition. The normal derivative of this function is: 


9G(x,y) 

9n 


9G(x,y) 

9r 


~ik H (2) 

2 "l , 


(2kr 




(30) 


which is bounded and continuous. 

For the interior region the elements used for the finite element analysis repre- 
sent a simple polar grid as shown in Figure 3. This necessitated two basic elements; 
a "triangular” element and a "rectangular” element. For this interior region the dis- 
placement was approximated by shape functions which are piecewise linear in both r 
and 0. This choice of shupe functions provided the necessary accuracy but kept the 
computations simple. 

The boundary elements of the halfspace for the functions <f> and \ are "line” 
elements. To preserve symmetry of the resulting matrices, these two functions used 
the same shape function. A piecewise constant shape function was chosen here, 
again because this provided the necessary accuracy while keeping computations simple. 
These elements were chosen to coincide with the boundary elements resulting from the 
interiors division. This is not necessary to the method but is more convenient 
computationally. 
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In the interior, the number of elements in the radial and angular directions was 
varied to provide a check of the convergence rate of the method. Although it 
resulted in elements of significantly different areas, the radial divisions were made 
uniform. This resulted in a natural placing of more elements where needed, near 
the center. The maximum number of elements used for this study was dictated by 
the limits of the computer used. 

The most computationally expensive parts of this method are the integrals 
involving the Green's function and its normal derivative. Because of their complexity, 
these integrals must be evaluated numerically. Gaussian quadrature was selected for 
this purpose. Because the integrand contains a logarithmic singularity, a modified 
Gaussian method, to take directly into effect this singularity, was used as presented 
by Harris and Evans [4]. In addition, care was taken to split the integral, at the 
singularity, into two integrals. The integration schemes used were tested on the 
natural logarithm function to determine the number of points necessary for accuracy. 
By using these techniques the amount of computation for each integral was kept to a 
minimum while returning the necessary accuracy. 

The symmetry of the full space problem was also used to reduce the number of 
computations. By constructing the matrices for the full space the actual number of 
integrals necessary to perform was reduced because of symmetry. Further symmetry 
considerations then allowed the matrices to be reduced to the halfspace problem. 

For the various matrix operations, such as solving equations and inverting 
matrices, the appropriate IMSL routines were used. 

IV. NUMERICAL RESULTS 


As stated previously, the primary objectives of this study were to determine 
the accuracy and convergence of the proposed formulation. Given the applications 
discussed in Section III a number of cases were considered. 

The first consideration was the number of elements needed to achieve a desired 
accuracy and the rate of convergence as the number of elements is increased. For 
this study, results were computed for five radial and five angular elements; 10 radial 
and 10 angular elements; and 20 radial and 20 angular elements. 


11 


The second consideration was the frequency of the free-field motion. As the 
frequency of the free-field motion increases, it was expected that the smoothness of 
the response would decrease. This indicates that more elements are needed to achieve 
a fixed accuracy for higher frequencies. The circular dimensionless frequencies con- 
sidered in this study were tt, 0.50tt, and 0.25ir. 

The third consideration was the angle of incidence of the free-field wave front. 
It was of interest to determine whether the angle of incidence has any effect on the 
solution. To do this, this study computed results for incidences of 0° and 60°. 

This represented wave fronts coming from directly below and from the left, respec- 
tively. 

The fourth consideration was the stiffness of the obstacle. The ability to vary 
the stiffness of the obstacle in relation to the stiffness of the halfspace is important 
to the usefulness of the formulation and any further extensions. For this study two 
cases were considered. An obstacle stiffer than the halfspace was studied using a 
ratio of mass densities (p"7p + ) equal to 1.50 and a ratio of shear wave velocities 

(g“/3 + ) equal to 2,00. A softer obstacle used a density ratio of 0.67 and a velocity 

ratio of 0.50. 

For all of the above applications and cases a closed form solution is available in 
terms of a series involving Bessel functions [2], All of these cases and their exact 
solutions were computed. The results for various points on the surface (9 = 0° and 
6 = 180°) were compiled into Tables 1 through 3. Both the real and imaginary parts 
of the solution are shown. For an incidence angle of 0° only those points where 
6 = 0° are shown since the response is symmetric. 


The relative errors for three surface points were computed and plotted against 
the number of elements for the various cases. This error is defined as: 


f 

error = 


exact - commuted 
exact 



The graphs of the errors are shown in Figures 4 through 9. 

As can be seen from the graphs, the convergence and accuracy of the formula- 
tion is very good. The rate of convergence is on the order of the square of the 
number of elements or better. The relative error itself is small even when only using 
20 elements. The effect of frequency generally follows the expected trend, that is, 
as the frequency increases more elements are needed to provide accurate answers. 
There is little effect on the accuracy of the formulation from the angle of incidence 
of the free-field motion or the relative stiffness of the obstacle and halfspace. 

Another concern was the viability of the formulation at or near the natural 
frequency of the obstacle when having the same material properties as the exterior 
region. At this natural frequency, the matrix [C] defined in equation (24) becomes 
ill-conditioned and difficult to invert. This is a particular problem with this study 
since inversion of the [C] matrix was required for the reduction done in equation 
(25). This inversion did, in fact, fail at the natural frequency of 0.7655ir. However, 
good results were obtained at a frequency of 0. 7639ir. Results for this frequency and 
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another close frequency are shown in Table 4. It is believed that this problem can 
be avoided by solving the full problem [equation (24)]. This could not be confirmed 
in this study, however, because of computer limitations. 


V. CONCLUSIONS 


The numerical results obtained indicate the formulation is very good. Accuracy 
is good even with few elements and the convergence rates are on the order of the 
number of elements squared or better. The formulation can be used for a wide range 
of incidence angles and frequencies. A variety of obstacle densities can be accommo- 
dated and the natural frequency of the halfspace can be closely approached. 

Now that the formulation has been numerically proven, a number of extensions 
for further study suggest themselves. An obstacle of inhomogeneous properties, a 
region of irregular shape, and a region of nonlinear properties are some possible 
extensions. This study provides a firm baseline for this further study. 
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Figure 5. 






























TABLE 1 



FREQUENCY * .25 n INCIDENCE ANGLE « 0° 



PjP+ = 

.67 

- ,50 


RADIUS 

5 ELEMENTS 

10 ELEMENTS 

20 ELEMENTS 

EXACT 

1.00 

1.1717 - 0.0106 

1.1320 - 0.0127 

1.1241 - 0.0129 

1.1218 - 0.0130 

0.80 

1.3671 - 0.0198 

1.3479 - 0.0220 

1.3444 - 0.0222 

1.3434 - 0.0222 

0.60 

1.5405 - 0.0278 

1.5338 - 0.0302 

1.5327 - 0.0304 

1.5324 - 0.0303 

0.40 

1.6787 - 0.0339 

1.6771 - 0.0365 

1.6769 - 0.0367 

1.6768 - 0.0366 

0.20 

1.7713 - 0.0379 

1.7685 - 0.0405 

1.7677 - 0.0407 

1.7674 - 0.0406 

0.00 

1.8209 - 0.0394 

1,8051 - 0.0419 

1.8003 - 0,0421 

1.7983 - 0.0420 


FREQUENCY = .50 n INCIDENCE ANGLE = 0° 



PJP+ " 

•67 

* .50 


RADIUS 

5 ELEMENTS 

10 ELEMENTS 

20 ELEMENTS 

EXACT 

1.00 

1.4287 - 0.1334 

1.2486 - 0.1654 

1.2172 - 0.1699 

1.2084 - 0.1710 

0.80 

1.0713 0.2401 

0.9629 0.2230 

0.9697 0.2214 

0.9661 0.2211 

0.60 

0.5390 0.7157 

0.5122 0.7145 

0.5099 0.7156 

0.5094 0.7161 

0.40 

- 0.0179 1.1896 

- 0.0109 1.1970 

- 0.0083 1.1988 

- 0.0074 1.1993 

0.20 

- 0.4351 1.5505 

- 0.4144 1.5520 

- 0.4105 1.5509 

- 0.4094 1.5503 

0.00 

- 0.5755 1.7429 

- 0.5650 1.7012 

- 0.5619 1.6856 

- 0.5605 1.6785 


FREQUENCY - 

= ,25 rr INCIDENCE ANGLE = 60° 



PjP + = 

■67 

= .50 


RADIUS 

5 ELEMENTS 

10 ELEMENTS 

20 ELEMENTS 

EXACT 

1.00 

0.4901 0.9179 

0.4675 0.9666 

0.4386 0.9894 

0.4010 1.0249 

0.80 

0.8362 0.8943 

0.8087 0.9093 

0.7939 0.9137 

0.7802 0.9179 

0.60 

' 1.1733 0.7559 

1.1533 0.7486 

1.1436 0.7454 

1.1336 0.7437 

0.40 

1.4662 0.5419 

1.4500 0.5268 

1,4424 0.5202 

1.4352 0.5143 

0.20 

1.6883 0.2828 

1.6732 0.2647 

1.6672 0.2553 

1.6623 0.2450 

0.00 

1.8274 0.0057 

1.8087 - 0.0171 

1 . 8021 - 0.0293 

1.7983 - 0.042 

- 0.20 

1.8405 - 0.2720 

1,8345 - 0.2982 

1.8331 - 0.3128 

1.8328 - 0.3272 

- 0.40 

1.7617 - 0.5325 

1.7614 - 0.5586 

1.7622 - 0.5743 

1.7635 - 0.5915 

- 0.60 

1.5915 - 0.7488 

1.5923 - 0.7778 

1.5939 - 0.7943 

1.5959 - 0.8129 

- 0.80 

1.3473 - 0.8902 

1.3435 - 0.9350 

1.3422 - 0.9557 

1.3428 - 0.9766 

- 1.00 

1.0580 - 0.9172 

1.0600 - 0.9884 

1.0455 - 1.0233 

1.0237 - 1.0714 


FREQUENCY = 

= .50 n INCIDENCE ANGLE = 60° 



u 

+ 

i 

,67 

= .50 


RADIUS 

5 ELEMENTS 

10 ELEMENTS 

20 ELEMENTS 

EXACT 

1.00 

- 0.8350 0.2944 

- 0.8821 0.2526 

- 0.9628 0.2343 

- 1.0670 0.2269 

0.80 

- 1.4689 1.4214 

- 1.5942 1.3720 

- 1.6468 1.3506 

- 1.6858 1.3345 

0.60 

- 1.7754 2.2891 

- 1.8995 2.2567 

- 1.9418 2.2458 

- 1.9703 2.2405 

0.40 

- 1.6925 2.6565 

- 1.7884 2.6518 

- 1.8177 2.6551 

- 1.8326 2.6631 

0.20 

- 1.2614 2.4187 

- 1.3040 2.4297 

- 1.3125 2.4424 

- 1.3115 2.4603 

0.00 

- 0.6277 1.6562 

- 0.5936 1.6530 

- 0.5765 1.6609 

- 0.5605 1.6785 

- 0.20 

0.1049 0.4709 

0.1593 0.5012 

0.1824 0.5196 

0.2013 0.5401 

- 0.40 

0.6459 - 0.6739 

0.7182 - 0.6506 

0.7436 - 0.6387 

0.7595 - 0.6270 

- 0.60 

0.8930 - 1.4961 

0.9498 - 1.4930 

0.9670 - 1.4922 

0.9724 - 1.4935 

- 0.80 

0.8273 - 1.7845 

0.8337 - 1.8073 

0 . 8253 — 1.8183 

0.8131 - 1.8324 

- 1.00 

0.5360 - 1.4826 

0.5303 - 1.5278 

0.4653 - 1.5522 

0.3734 - 1.5880 


TABLE 2 



FREQUENCY * 

,25* 

INCIDENCE ANGLE * 

0° 




PjP + - 1.50 


= 2.00 




RADIUS 

5 ELEMENTS 

10 ELEMENTS 

20 ELEMENTS 

EXACT 

1,00 

0.8880 

0.1357 

0.8863 

0.1339 

0.8858 

0.1336 

0.8856 

0.1336 

0.80 

0.8840 

0.1380 

0.8837 

0.1361 

0.8835 

0.1358 

0.8835 

0.1358 

0.60 

0.8816 

0.1397 

0.8819 

0.1378 

0.8819 

0.1375 

0.8819 

0.1375 

0.40 

0.8803 

0.1410 

0.8807 

0.1391 

0.8808 

0.1388 

0,8808 

0.1368 

0.20 

0.8799 

0.1418 

0.8801 

0.1398 

0.8802 

0.1395 

0.8803 

0.1395 

0.00 

0.8807 

0.1421 

0.8802 

0.1401 

0.8800 

0.1398 

0.8800 

0.1398 


FREQUENCY = 

.50* 

INCIDENCE ANGLE - 

0° 




PjP+ = 1-B0 


- 2,00 




RADIUS 

5 ELEMENTS 

io elements 

20 ELEMENTS 

EXACT 


1.00 

0.5750 

- 0.0041 

0.5729 

- 0 . GU 52 

0.5714 

- 0.0052 

0.5706 

- 0.0051 

0.80 

0.5549 

0.0065 

0.5564 

0.0046 

0.5563 

0.0044 

0.5561 

0.0044 

0.60 

0.5421 

0.0148 

0.5451 

0.0124 

0.5456 

0.0121 

0.5456 

0.0120 

0.40 

0.5346 

0.0206 

0.5379 

0.0181 

0.5385 

0.0177 

0.5386 

0.0176 

0.20 

0.5313 

0.0240 

0.5341 

0.0215 

0.5345 

0.0211 

0.5346 

0.0210 

0.00 

0.5327 

0.0249 

0.5335 

0.0225 

0.5334 

0.0222 

0.5333 

0.0222 


FREQUENCY = 

,25* 

INCIDENCE ANGLE = 

60° 




pjp* " 1.50 


= 2.00 




RADIUS 

5 ELEMENTS 

10 ELEMENTS 

20 ELEMENTS 

EXACT 


1.00 

0.8750 

0.3010 

0.8728 

0.2936 

0.8711 

0.2907 

0.8683 

0.2894 

0.80 

0.8828 

0.2775 

0.8809 

0.2694 

0.8800 

0.2652 

0.8793 

0.2610 

0.60 

0.8876 

0.2510 

0.8866 

0.2418 

0.8862 

0.2368 

0.8860 

0.2315 

0.40 

0.8888 

0.2228 

0.8884 

0.2127 

0.8883 

0.2073 

0.8883 

0.2016 

0.20 

0.8861 

0.1936 

0 . 8861 

0.1828 

0.8861 

0.1770 

0.8862 

0.1637 

0.00 

0.8798 

0.1644 

0.8797 

0.1525 

0.8798 

0 . 1461 

0.8800 

0.1398 

- 0.20 

0.8691 

0.1346 

0.8690 

0.1217 

0.8691 

0.1149 

0.8693 

0.1155 

- 0.40 

0.8543 

0.1039 

0.8542 

0.0905 

0.8544 

0.0833 

0.8545 

0.0763 

- 0.60 

0.8357 

0.0733 

0.8355 

0.0593 

0.8355 

0.0517 

0.8355 

0.0445 

- 0.80 

0.8136 

0.0436 

0.8131 

0.0286 

0.8127 

0.0205 

0.8123 

0.0122 

- 1.00 

0.7889 

0.0159 

0.7885 

0.0006 

0.7875 

*“ 0.0087 

0.7852 

- 0.0197 


FREQUENCY = 

,50* 

INCIDENCE ANGLE = 

60° 




P //> = 1,50 

“* T 


= 2,00 




RADIUS 

5 ELEMENTS 

10 ELEMENTS 

20 ELEMENTS 

EXACT 


1.00 

0.5944 

0.2505 

0.5941 

0.2412 

0.5930 

0.2374 

0.5883 

0.2345 

0.80 

0.5973 

0.2141 

0.5991 

0.2052 

0.6009 

0.2006 

0.6032 

0.1960 

0.60 

0.5914 

0.1740 

0.5969 

0.1648 

0.6006 

0.1597 

0.6047 

0.1544 

0.40 

0.5753 

0.1318 

0.5834 

0.1221 

0.5881 

0.1167 

0.5932 

0.1114 

0.20 

0.5486 

0.0886 

0.5582 

0.0784 

0.5635 

0.0726 

0.5692 

0.0676 

0.00 

0.5126 

0.0462 

0.5221 

0.0347 

0.5275 

0.0285 

0.5333 

0.0222 

- 0.20 

0.4659 

0.0038 

0.4751 

- 0.0085 

0.4805 

- 0.0150 

0.4862 

- 0.0225 

- 0 . 40 

0.4095 

- 0.0380 

0.4183 

- 0.0505 

0.4234 

- 0.0571 

0.4288 

- 0.0641 

- 0.60 

0.3450 

- 0.0777 

0.3528 

- 0.0901 

0.3573 

- 0.0968 

0.3619 

- 0.1035 

- 0.80 

0.2744 

- 0.1140 

0.2803 

- 0.1264 

0.2834 

- 0.1332 

0.2865 

- 0.1403 

- 1.00 

0.2006 

- 0.1454 

0.2063 

- 0.1572 

0.2071 

- 0.1645 

0.2035 

- 0.1728 
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TABLE 3 


FREQUENCY = 1,00* INCIDENCE ANGLE = 0° 
pJP A * .67 * .50 


RADIUS 

5 ELEMENTS 

10 ELEMENTS 

20 ELEMENTS 

EXACT 

1.00 

2.6999 0.4648 

1.5254 - 0.0488 

1.4244 - 0.0224 

1.4034 - 0.0037 

0.80 

1.6630 0.9319 

1.1053 0.4798 

1.1905 0.5251 

1.2366 0.5571 

0.60 

- 0.1468 1.0585 

- 0.1258 0.7668 

0.0042 0.7883 

0.0607 0.8082 

0.40 

- 1.4253 0.6158 

- 1.2088 0.4260 

- 1.1236 0.4102 

- 1.0900 0.4099 

0.20 

- 1.8934 - 0.1374 

- 1.6905 - 0.2832 

- 1.6331 - 0.3128 

- 1.6127 - 0.3216 

0.00 

- 2.3700 - 0.5224 

- 1.8809 - 0.6694 

- 1.7531 - 0.6850 

- 1,7050 - 0.6859 


FREQUENCY - 

- 1.00* INCIDENCE ANGLE = 60° 



PJP+ = 

.67 

« .50 


RADIUS 

5 ELEMENTS 

10 ELEMENTS 

20 ELEMENTS 

EXACT 

1.00 

- 1.1846 0.8262 

- 0.5685 0.4656 

- 0.6292 0.3569 

- 0.6971 0.1962 

0.80 

0.2819 0.5360 

0.6327 0.5367 

0.6107 0.6343 

0.6021 0.7031 

0.60 

0.8912 - 0.4740 

0.9871 - 0.3300 

0.9713 - 0.2135 

0 . S 546 - 0.1176 

0.40 

- 0.3372 - 1.4100 

- 0.3591 - 1.5494 

- 0.3842 - 1.5829 

- 0.4230 - 1.5687 

0.20 

- 2.1465 - 1.4241 

- 1.9791 - 1.7901 

- 1.9613 - 1.9122 

- 1.9885 - 1.9634 

0.00 

- 2.3429 - 0.4825 

- 1,8655 - 0.6399 

- 1.7454 - 0.6691 

- 1.7050 - 0.6859 

- 0.20 

- 0.1101 0.6409 

0.1856 0.8438 

0.2620 0.9304 

0.3148 0.9629 

- 0.40 

1.9918 0.7143 

1.8290 1.1634 

1.7830 1.2555 

1.7968 1.2719 

- 0 . 60 

2.0383 - 0.0450 

1.3603 0.2853 

1.1410 0.2157 

1.0631 0.1366 

- 0.80 

0.0976 - 0.8284 

- 0.5279 - 0.6316 

- 0.7781 - 0.7939 

- 0.8945 - 0.9124 

- 1.00 

- 1.7652 - 0.9353 

- 1.4990 - 0.7225 

- 1.5759 - 0.7304 

- 1.6672 - 0.6452 


FREQUENCY * 

= 1,00* INCIDENCE ANGLE = 0° 



PjP + = 

1.50 J 3 JJ 3 + 

- 2.00 


RADIUS 

5 ELEMENTS 

10 ELEMENTS 

20 ELEMENTS 

EXACT 

1.00 

0.2325 - 0.3428 

0.2463 - 0.3417 

0.2431 - 0.3396 

0.2340 - 0.3387 

0.80 

0.1802 - 0.3015 

0.1897 - 0.3047 

0.1908 - 0.3045 

0.1890 - 0.3044 

0.60 

0.1525 - 0.2679 

0.1596 - 0.2729 

0.1612 - 0.2736 

0.1614 - 0.2738 

0.40 

0.1390 - 0.2440 

0.1445 - 0.2486 

0.1458 - 0.2496 

0.1461 - 0.2498 

0.20 

0.1333 - 0.2307 

0.1379 - 0.2338 

0.1387 - 0.2344 

0.1389 - 0.2345 

0.00 

0.1343 - 0.2314 

0.1367 - 0.2301 

0.1369 - 0.2295 

0.1368 - 0.2293 


FREQUENCY * 

= 1,00* INCIDENCE ANGLE = 60° 



PjP + ~ ' 

1,50 ftjft* 

= 2 . 1)0 


RADIUS 

5 ELEMENTS 

10 ELEMENTS 

20 ELEMENTS 

EXACT 

1.00 

0.5797 0.0607 

0.5966 0.0309 

0.6050 0.0153 

0.6036 - 0.0107 

0.80 

0.5697 - 0.0346 

0.5835 - 0.0500 

0.5929 - 0.0604 

0,6013 - 0.0729 

0.60 

0.5107 - 0.1136 

0.5260 - 0.1194 

0.5352 - 0.1257 

0.5448 - 0.1340 

0.40 

0.4107 - 0.1735 

0.4253 - 0.1752 

0.4333 - 0.1789 

0.4416 - 0.1843 

0.20 

0 . 2794 - 0.2126 

0.2902 - 0.2126 

0.2959 - 0.2143 

0.3016 - 0.2172 

0.00 

0.1292 - 0.2299 

0.1324 - 0.2288 

0.1345 - 0.2288 

0.1368 - 0.2293 

- 0.20 

- 0.0319 - 0.2252 

- 0.0368 - 0.2225 

- 0.0387 - 0.2210 

- 0.0401 - 0.2192 

- 0.40 

- 0.1937 - 0.1995 

- 0.2058 - 0.1951 

- 0,2115 - 0.1917 

- 0.2167 - 0.1874 

- 0.60 

- 0.3441 - 0.1555 

- 0.3634 - 0.1478 

- 0.3726 - 0.1424 

- 0.3819 - 0.1355 

- 0.80 

- 0.4743 - 0.0975 

- 0.4994 - 0.0835 

- 0.5119 - 0.0751 

- 0.5265 - 0.0646 

- 1.00 

- 0.5787 - 0.0319 

- 0.6017 - 0.0096 

- 0.6170 0.0017 

- 0.6438 0.0248 
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TABLE 4 


RADIUS 

1.00 

0.80 

0.60 

0.40 

0.20 

0.00 


R ADIUS 

1.00 

0.80 

0.60 

0.40 

0.20 

0.00 

- 0.20 

- 0.40 

- 0.60 

- 0.80 

- 1.00 


RADIUS 

1.00 

0.80 

0.60 

0.40 

0.20 

0.00 


RADIUS 

1.00 

0.80 

0.60 

0.40 

0.20 

0.00 

- 0.20 

- 0.40 

- 0.60 

- 0.80 

- 1.00 


FREQUENCY * 7639* INCIDENCE ANGLE * 0° 

p J p + 85 -67 s *50 

5 ELEMENTS 10 ELEMENTS 20 ELEMENTS EXACT 


1.2885 0.7557 

0.6008 

0.5329 

0.4699 0.4502 

0.4322 

0.4191 

1.7459 2.1061 

0.8836 

2.2977 

0.6170 2.2960 

0.5240 

2.2857 

1.5807 2.6912 

0.7885 

3.1921 

0.4843 3.2716 

0.3731 

3.2874 

0.8648 2.4518 

0.3467 

2.9921 

0.1232 3.1032 

0.0393 

3.1337 

0.0481 1.8992 

- 0.1567 

2.2912 

- 0.2559 2.3843 

- 0.2947 

2.4139 

- 0.2416 1.9265 

- 0.3773 

2.0175 

- 0.4188 2.0315 

- 0.4336 

2.0327 

FREQUENCY * 

= .7639* 

INCIDENCE ANGLE = 60° 



PjP+ = 

.67 


= .50 



5 ELEMENTS 

10 ELEMENTS 

20 ELEMENTS 

EXACT 

0.1824 - 0.0520 

0,5011 

- 0.1965 

0.5007 - 0.2320 

0.4419 

- 0.3117 

- 0.7646 - 1.7509 

- 0.2770 

- 2,1918 

- 0.0889 - 2.3035 

0.0121 

- 2.3884 

- 1.5071 - 2.4883 

- 1.0001 

- 2.9706 

- 0.7766 - 3.1009 

- 0.6631 

- 3.1864 

- 1.6157 - 1.7282 

- 1.3282 

- 1.9749 

- 1,1970 - 2.0368 

- 1,1350 

- 2.0708 

- 1.0528 0.0848 

- 1.0578 

0.1435 

- 1.0509 0.1761 

- 1.0525 

0.2039 

- 0.2289 1.8186 

- 0.3652 

1.9493 

- 0.4119 1.9956 

- 0.4336 

2.0327 

0.4553 2.2009 

0.3285 

2.1436 

0.3190 2.1440 

0.3292 

2.1542 

0.4445 1.1071 

0.5636 

0.7313 

0.6704 0.6352 

0.7445 

0.5848 

- 0.0427 - 0.6358 

0.2801 

- 1.1881 

0.4672 - 1.3377 

0.5750 

- 1.4226 

- 0.5517 - 1.8617 

- 0.2191 

- 2.2071 

- 0.0723 - 2.2967 

0.0011 

- 2.3533 

- 0.6603 - 1.8214 

- 0.4074 

- 1.6962 

- 0.4329 - 1.6356 

- 0.5248 

- 1.5736 

FREQUENCY = 

- .7480* 

INCIDENCE ANGLE = 0° 



PJP+ " 

.67 


= .50 



5 ELEMENTS 

10 ELEMENTS 

20 ELEMENTS 

EXACT 

1.4014 0.7395 

0.7765 

0.6111 

0.6452 0.5666 

0.6046 

0.5495 

1.8342 1.8155 

1.2719 

2.0665 

1.0756 2.1275 

1.0040 

2.1448 

1.6919 2.2815 

1.2144 

2.7801 

1.0143 2.9148 

0.9374 

2.9573 

0.9262 2.1212 

0.6445 

2.6158 

0.5058 2.7528 

0.4511 

2.7974 

0.0697 1.7310 

- 0.0322 

2.0636 

- 0.0911 2.1530 

- 0.1149 

2.1824 

- 0.2267 1.7944 

- 0.3242 

1.8570 

- 0.3524 1.8639 

- 0.3617 

1.8621 

FREQUENCY = 

= .7480* 

INCIDENCE ANGLE = 60° 

f 


Pjp + = 

,67 


= .50 



5 ELEMENTS 

10 ELEMENTS 

20 ELEMENTS 

EXACT 

0.1479 - 0.0826 

0.4103 

- 0.2591 

0.4024 - 0.3148 

0.3392 

- 0.4012 

- 0.8543 - 1.6355 

- 0.5289 

- 2.0535 

- 0.3988 - 2.1883 

- 0.3282 

- 2.2870 

- 1.6098 - 2.3271 

- 1.2809 

- 2.7524 

- 1.1303 - 2.8988 

- 1.0544 

- 2.9958 

- 1.6969 - 1.6638 

- 1.5085 

- 1.8667 

- 1.4222 - 1.9385 

- 1.3828 

- 1.9807 

- 1.0934 - 0.0005 

- 1.0945 

0.0607 

- 1.0901 0.0850 

- 1.0925 

0.1067 

- 0.2196 1.6905 

- 0.3166 

1.7924 

- 0.3480 1.8299 

- 0.3617 

1.8621 

0.5098 2.2102 

0.3617 

2.1399 

0.3431 2.1275 

0.3491 

2.1296 

0.5258 1.3568 

0.5114 

0.9991 

0.5652 0.8854 

0.6136 

0.8226 

0.0252 - 0.2368 

0.1475 

- 0.7638 

0.2579 - 0.9382 

0.3291 

- 1.0405 

- 0.5396 - 1.5412 

- 0.3682 

- 1.9033 

- 0.2769 - 2.0252 

- 0.2283 

- 2.1019 

- 0.7173 - 1.7779 

- 0.5017 

- 1.7239 

- 0.5288 - 1.6946 

- 0.6191 

- 1.6529 
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APPENDIX A 


Substituting the shape functions [equation (23)] into equation (22) is done 
after grouping certain terms. 

1st Term 


~ I {y~(w <5 w + w <$w ) - w 2 p~ w <$ w} dx 

•'SI X 1 X 1 x 2 x 2 

= -/ | { 6 w } T [N ] [N ] T lW} + {<S w> T [N ] jAN 1 T {w> 
fi ( X 1 X 1 x 2 x 2 


w 2 {6 w} T [N] p"[N] T {w}> dx 


Let [K] = -/{[ N v ] y-[N ] T + [N v ] ii~[N v ] T - <o 2 [N] T p“[N] } dx 
Q X 1 X 1 x 2 x 2 


then the 1st term becomes 


,T 


(5W> A [K] {w} . 


2nd Term 


/ w"i 


1/2 f w~(x) (Sp“ w "(x) dS 
r 


/ 


1/2 J { 6 A } T [Q] [N r ) T {w p } dS 


/ 


Let [A] = 1/2 J [Q] [NJ T dS 


then the 2nd term becomes 




V’J 


3rd Term 


r ( + + f 3G<x,y) > 

1/2 J 1 1/2 pi w“(x) + u J — w"(jr) dS^ J H(x) dS 


= 1/2 jf j 1/2 / {6<j>} T [Q] [N r 3 T {Wp } dS 


+ p 


f m 8G(X,jr) rn ) 

y CQi tNp) 1 (w r }ds z ds 




Let [B] = 1/2 y + 1 1/2 J [Q] [N r l T dS + 


•ffw pT§r] ,M i- ,T 


dS y dS 


then the 3rd term becomes 


in) [B] {Wp> . 


4th Term 


- 1/2 I J G(x,y) p“ w “(y) 6<j> dSy dS 
r r 


/ / G (£> 

r r 

a 

a 

r r 

then the 4th term becomes 


= 1/2 J J (6 0 [Q] [G(x,y>] [Q] x {A } dSy dS 
r r 


Let [C] = - 1/2 J J [Q] [G(x,y>] [Q} x dSy dS 

r r 


{«*} [C] {X} 


5th Term 


/ 


/ 


]i~ w 6 w dS - 1/2 J p“ w _ 6 w" dS = J {6w„ } T [N _ ] [Q] T {X} dS 

p ** p I* p ■* 


/ 


/ 


- 1/2 j {fiw } T [N ] [Q] T {X} dS 

r 


k 

>] 
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then the 5th term becomes 


{Sw r ) T CA] T {a} . 


6th Term 




■f 


r 3G(x,£) 

1/2 / \U2 \l' <j)(x) + J — ~ - 


4><X> ds X 




= 1/2 J {1/2 / {6w p } + [N r ] [Q] T 4} 


+ P + / 4w r ) T [N p ] 


3G(x,y) 


Let [B] T = 1/2 y + 


1 1/2 1 


3 nx 


,T 


[N r ] [Q] A dS + / [N ] 

r r J v r r L 


[Q] T 4) dSy 

/ IN r ] f 


then the 6th term becomes 


{<Sw r } T [B] T . 


7th Term 


// 


“1/2 J J G(x,y) <j>(y) 5 y“ w “(x) dSy dS 

r r 


- 1/2 


//tm T 


[Q] [G(x,y>] [Q] A 4) dSy dS 


r r 


Let [C] T = - 1/2 If [Q] [G(x,y>] [Q] T dSy dS 

r r 


then the 4th term becomes 


<5vv - dS 

| dS 

3G(x,y) 
3 nx 


8th Term 


/ 


- 1/2 J w°(x) fiy~ w n “(x) dS ~ - 1/2 


/ { 6 X } T 


[Q] w° dS 


Let {ny - 1/2 / [Q] w° dS 

r 

then the 8th term becomes 


- {6A} 1 {fX } . 


9th Term 


1/2 - / / w 0 6 w“ dS = 1/2 y + / {6 w„ } T [N P J w ° dS 


Let {f r > = - 1/2 y Jf [N f ] w n ° dS 


/ 


then the 9th term becomes 


{<Sw p } A {f r > . 


10th Term 




1/2 y J w° 6(J> dS = - 1/2 y ' J {6<J.} A [Q] w° dS 

r *r 


f {6 d. > T 


Let {f^} = y 1/2 J [Q] w° dS 


/ 


then the 10th term becomes 


{f + } . 
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Equation (22) now becomes: 

{6 w} T ([KJ {w}) + {6w r ) T ([B] (*) + [A] T {X} - {f r » 

+ {6 X } T ([A] <Wp } + [C] T {<#>> - {fX}) + (6 <#>> T ([B] {Wp} 

+ [C] {X} - {f<j>}) = 0 * 

For arbitrary (<5w>, {6 w }, { 6 X > , and {6^}; quantities multiplying them must each 
be zero . 
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APPENDIX B. 


PROGRAM LISTING 


C INTRODUCTION 
C ~ 

c 

c This program solves, using finite element techniques, a soil structure 
C interaction problem defined in polar coordinates. The structure is a 
C semicylindericai deposit of soil at the surface of an elastic halfspace, 
c The properties of this deposit can be diferent from those of the 
C half space and can also represent a linear change of shear modulus with 
C depth. The forcing function of this problem is an incoming antiplane 
c shear wave. Th computed results are the steady state antiplane 
C displacements. The interior deposit is represented by finite element 
C techniques while the exterior half space is represented at the boundary 
C of the deposit bye the use of Green's functions. Thus only finite 
C element storage requirements are for the deposite and no additional 
C finite element manipulations are necessary. 

C 

C PROGRAM STRUCTURE 

C — 

c 

c The program is stuctUred with a main driving program with various 
C levels of subroutines. A schematic of this structure is shown here. 

C 


c 

LEVEL 0 

LEVEL 1 LEVEL 2 

LEVEL 3 

LEVEL 4 






c 





c 


|— MESH 



c 


I 1 — RECT 



c 


| ASSM i 



c 


I ! — TRI 



c 


! 



c 


1 j — GAMMA 



c 


|— - FORCE ! 



c 


| |— PHI LAM 



c 


1 



c 

MAIN 

| I — BMAT 



c 


1 1 
1 1 

NORM — 


c 


I i~ AMAT — 

. 

DERIV 

c 


1 1 
1 1 

NORM2 — 

■ 

c 


!--- HSPACE — !— ATEST 



c 


I ! 

DIAG — 

. . : • . 

c 


I I— CMAT 


FUNC 

c 


! ! 

DIAG2 — 

- 

' 

c 


1 !-- CTEST 



c 


j KTEST 




C 

C A description of what each subroutine does is presented at the 
C beginning of that subroutine. 

C 

C INPUT DATA 

C — 

c 

C The input data for this program is contained in two data files. 

C F0R24.DAT contains the Gauss integration information and should not be 
c changed without a series of tests to check convergence. F0R25.DAT 
c contains all the other information necessary for the operation of this 
C program. The unit numbers for these files are 24 and 25 respectively. 

C An example of the data contained in this file is shown here: 

C 

c 

c AMAT | 

PRECEDING PAGE BEAN1C NOT FfBMED 31 

*«Ll2l2.«i£(BtfKAUif Dunk 




C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


CMAT 

HSPACE 

MESH 

START — 
UNIFORM 


| — debug switches (subroutine names) 


signals end of debug switches (must always be used) 
radial division method (UNIFORM or GRADIENT) 


5.76 0.0 0.5 0.5 1.0 1.0 0.5 

10 10 

3.14159 1.0 1 

1.0 1.0 1.0 1.0 1.0 


-frequency squared 
incidence angle (radians) 
shear modulus at depth 
shear modulus at surface 
density of deposit 
shear modulus of halfspace 
jgeometric shape factor ( .5) 

-number of radial divisions 
_number of angular divisions 

-angular sweep (radians) 
radial range 

-number of boundary elements 
per angular element 

boudary node displacements 
if KTEST is to be run. 


In use this file would look like this? 


AMAT 

CMAT 

HSPACE 

MESH 

START 

UNIFORM 

5.76 0.0 0.5 0.5 1.0 1.0 0.5 
10 10 

3.14159 1.0 1 

1.0 1.0 1.0 1.0 1.0 


OUTPUT DATA 


All output is directed to the unit number 26. 


MAIN PROGRAM 


The main or driving program has four major duties. The first is to 
read in the data which must be supplied by the user. A more detailed 
description of this data is presented in the "Input Data" section. The 
second duty of this subprogram is to set the bollion values for the 
various debug and method options available in this program. The 
third duty is the definition of matrix sizes. Each matrix is defined 
here using an actual number which represents the maximum size 
allowed. This limits the size of structure whcih can be analysed. 

Using the data which is read in the actual sizes of the matrices are 
calculated. These values are then used throughout the program, 
computed matrix sizes are passed in COMMON to all subroutines while 
the matrices themselves are passed as needed through the arguments of 
the various subroutines. An expansion of the size of structure the 
program can handle can be done 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

0 

c 

c 

c 


c 

c 

c 

G 

c 

G 

c 


from this subroutine alone since this is the only place where these 
matrices are defined using actual numbers. 

However, if the size of the individual elements used is changed then 
the subroutines ASSM, RECT, and TRI must also be changed and their 
accompanying matrices. 

The final duty of this subprogram is to control the execution of 
the entire program. This is done by calling vario;;." r’ib r outines which 
do specific jobs such as constructing particular matrices or solving 
the equations. 


Declare in COMMON all variables required by more 

than one C subroutine. These variables are divided into four 

catagoriess 

1. SIZE This common block contains variables which define 

array dimensions. 

2. PROB This common block contains variables which define 

problem parameters. 

3. MAP This common block contains variables which define 

the finite element mesh. 

4. dbug This common block contains variables which define 

which debug switches are desired. 

Arrays are not passed through COMMON but passed as arguments to the 
subroutines. 


COMMON /SIZE/ NUMANG , NUMRAD , NUMELM , INTNOD , NUMNOD ,CODIAG , 

1 NUMEQN , LENGTH ,WIDTH, SPACE , BNDELM , BCOLS , ALENG , AWIDE, ATMPL 

INTEGER NUMANG , NUMRAD , NUMELM , INTNOD , NUMNOD , CODI AG , NUMEQN , 

1 LENGTH .WIDTH, SPACE , BNDELM , BCOLS , ALENG , AW I DE , ATMPL 

COMMON /PROB/ W2 , GA , GB , RO , EXTG , ALPHA , THETAO 
REAL W2 , GA , GB , RO , EXTG .ALPHA, THETAO 

COMMON /MAP/ SWEEP, RANGE, DRAD, DANG, BDANG1 , BDANG2 , RATIO, GRADNT 
INTEGER RATIO 

REAL SWEEP , RANGE , DRAD , DANG , BDANG1 , BDANG2 
LOGICAL GRADNT 

common /dbug/ dmain,dmesh,dassm,dktest,dforce, 

1 dhspac , drect , dtri , dgamma , df ilam, damat , dbmat , 

1 dcmat,dnorm,ddiag,ddiag2,dfunc,dderiv, detest, 

1 datest 

logical dmain , dmesh , dassm , dktest , df orce , dhspac , 

1 drect , dtri , dgamma , df ilam , damat , dbmat , demat , dnorm , 

1 ddiag,ddiag2,dfunc,dderiv, detest, datest 


Declare and dimension arrays necessary. Only those arrays which are 
dependent on structure size are dimensioned here. This is the only 
location where actual numbers are used for these arrays. An increase 
of structure size can be handled from here alone. Arrays which are 
dependent on element size are dimensioned in RECT, TRI, and ASSM. 

INTEGER GEOM (400,4) 

REAL PROP (400, 5) ,K(421,45) ,0UTPUT(421) ,KOG(400,2l) , 

1 B (400 , 21) ,KOO(400,45) ,XL(400,23) ,WG(2l) ,WO(400) , 

2 EXACT (400) , BMATRX (20 , 2) 

COMPLEX FG(21) ,FPHI (21) ,C(20,20) ,CTEMP(20) ,FD(21) ,W(21) , 

1 CINVER(20,20) .ATRANS (21,20) ,ATEMP(2,21) ,ATEMP2 (80,40) , 

1 ATEMP3 (80,40) ,D(21,21) ,TEMPD{21,21) ,WHOLE{21,21) , 

1 LAM (20) 

Declare variables required in this main program. 

LOGICAL OK 
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INTEGER debug, METHOD 
REAL FREQ,G 
C 

C Define whether debug is required or not. 
C 

C Initialize, 
c 


dmain * . 

false. 

dmesh - . 

false. 

dassm = . 

false. 

dktest - 

.false. 

dforce * 

.false. 

dhspac = 

.false. 

dbmat - 

.false. 

drect = 

.false. 

dtri = 

.false. 

dgamma = 

.false. 

dfilam = 

.false. 

damat = 

.false. 

dbmat = 

.false. 

dcmat - 

.false. 

dnorm = 

.false. 

ddiag = 

•false. 

ddiag2 = 

.false. 

dfunc = 

.false. 

dderiv = 

.false. 

detest = 

.false. 

datest = 

•false. 


Read in debugs desired. 

5 read (25, 7) debug 

7 format (a5) 

Echo debugs desired. 

write (5, 8) debug 

8 format (2X,a5) 

C 

C Set switches for debugs desired. 'START' indicates no more 
c switches will be requested. 

C 

if (debug. eq. 'MAIN' )dmain = .true, 
if (debug. eq. 'MESH' )dmesh = .true, 
if (debug. eq. 'ASSM’ )dassm = .true, 
if (debug. eq. 'KTEST' )dktest = .true, 
if (debug. eq. 'FORCE Odforce - .true, 
if (debug. eq. 'HSPAC' )dhspac = .true, 
if (debug. eq. 'RECT' )drect = .true, 
if (debug. eq. 'TRI ’ )dtri * .true, 
if (debug. eq. 'GAMMA' )dgamma = .true, 
if (debug, eq. 'PHILA' jdfilam 55 .true, 
if (debug. eq. 'AMAT' )damat - .true, 
if (debug. eq. 'BMAT' ) dbmat = .true, 
if (debug. eq. 'CMAT' jdcmat = .true, 
if (debug. eq. ' NORM' )dnorm = .true, 
if (debug. eq. 'DIAG'Jddiag = .true, 
if (debug. eq. 'DIAG2')ddiag2 = .true, 
if (debug. eq. 'FUNC' )dfunc = .true, 
if (debug. eq. 'DERIV' )dderiv » .true. 
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if (debug. eg. 'CTEST') detest * .true, 
if (debug. eq. 'ATEST' jdatest * .true, 
if (debug. eq. ’START') goto 10 
goto 5 

10 continue 

Define desired radial division method. If GRADIENT is desired radial 
divisions are chosen to balance element length and width. If UNIFORM 
is requested radial divisions are equal segments. 

Initialize. 

GRADNT * .FALSE. 

Read in desired method. 

READ (2 5, 11) METHOD 

11 FORMAT (1A5) 

Echo method desired. 

WRITE (5 ,12) METHOD 

12 FORMAT (2X,1A5) 

Set switch. 

IF (METHOD. EQ. ' GRADI ' ) GRADNT = .TRUE. 


Read in incident wave information and soil parameters. 
READ (25 , * ) W2 , THETAO , G , GB , RO , EXTG , ALPHA 


Compute frequency from input squared frequency. 

FREQ = SQRT (W2) 

Read in desired number of angular and radial divisions. 

READ (25, * ) NUMANG, NUMRAD 

Read in angle sweep of the structure (o-?radians) and the radius 
range of the structure (0-?) . 

READ (25 , « ) SWEEP , RANGE, RATIO 

Calculate the rate of change of soil stiffness. 

GA = (G-G3) /RANGE 

Calculate array dimensions and make checks for program capacity. 
OK * .TRUE. 

C# IF (RATIO. GE. 1. AND. RATIO. LE. 5) GOTO 15 

C# WRITE(5,800) 

C# OK = .FALSE. 

15 BNDELM = NUMANG/RATIO 
C# I = MOD (NUMANG, RATIO) 

C# I F ( I . EQ . 0 . AND . BNDELM . LE . 15 ) iL’OTO 16 

C# WRITE (5, 150) 

C# OK = .FALSE. 


16 

BCOLS * RATIO + 1 

C# 

IF(BC0LS.LE.6) GOTO 17 

c# 

WRITE(5,150) 

c# 

OK * .FALSE. 

17 

AW IDE « BNDELM*2 

C# 

IF (AWIDE.LE.30) GOTO 18 

c# 

WRITE(5,150) 

c# 

OK « .FALSE. 

18 

ALENG * BNDELM*RATIO*4 

C# 

IF (ALENG . LE . 300) GOTO 19 

c# 

WRITE(5,150) 

c# 

OK " .FALSE. 

19 

ATMPL ■ RATIO* 2 

c# 

IF(ATMPL.LE.10)GOTO 20 

c# 

WRITE(5,150) 

c# 

OK * .FALSE. 

20 

NUMELM * NUMANG*NUMRAD 

c# 

IF (NUMELM . LE. 225) GOTO 21 

C# 

WRITE(5,100) 

c# 

OK « .FALSE. 

21 

INTNOD * NUMANG+1 

c# 

IF ( INTNOD . LE . 16) GOTO 30 

c# 

WRITE(5,200) 

c# 

OK = .FALSE. 

30 

NUMNOD = (INTNOD) *NUMRAD+1 

c# 

1F(NUMN0D.LE.241) GOTO 40 

c# 

WRITE(5, 300) 

c# 

OK = .FALSE. 

40 

CODIAG ~ NUMANG + 2 

c# 

IF(C0DIAG.LE.17) GOTO 50 

c# 

WRITE(5,400) 

c# 

OK = .FALSE. 

50 

NUMEQN = NUMNOD-INTNOD 

c# 

IF (NUMEQN . LE . 225) GOTO 60 

c# 

WRITE(5,500) 

c# 

OK = .FALSE. 

60 

LENGTH = NUMNOD* (NUMNOD+l), 

c# 

IF ( LENGTH. LE. 29161) GOTO 7< 

c# 

WRITE(5,600) 

c# 

OK “ .FALSE. 

70 

WIDTH = 2«CODIAG+l 

C 

SPACE = CODIAG+1 


c If checks for program capacity are allright continue otherwise abort 
C program execution. 

C 

IF (OK) GOTO 80 
WRITE(5,700) 

STOP 

C 

C Generate finite element mesh. 

C 

80 CALL MESH ( PROP, GEOM) 

C 

C Assemble structure stiffness matrix, 

c 

CALL ASSM( PROP, GEOM, K, OUTPUT) 

C 

C Determine applied forces. 

C 
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CALL F0RCE(FG*FPHI) 


Generate matrix defining half space. 

CALL HSPACE ( ATRANS , ATEMP , ATEMP2 , ATEMP3 , BMATRX , 

1 C , CTEMP , FG , FPH I , FD , OUTPUT , CIN VER , TEMPO , D) 

I test of K matrics is desired call KTESI subroutine. 

IF ( DKTEST) CALL KTEST ( K , WO , WG , KOG , B * KOO , XL , OUTPUT , EXACT) 

Solve for displacements. 

CALL SOLVE (K , D , FG , FD , KOG , KOO , B , WHOLE , SLTN , W , XL , OUTPUT , 

1 LAM , TEMPD , Cl N VER , FPH I ) 

output computed values. 

CALL OUT(SLTN,W,FG,FD,B,LAM) 

Output information supplied to program by user. 

WRITE (26 , 1000) SWEEP , NUMAHG , RANGE , NUMRAD , ALPHA , RATIO 

IF(GRADNT) WRITE (26, 1500) 

IF(. NOT. GRADNT)WRITE{26, 1600) 

WRITE { 26 , 2 000 ) W2 , FREQ , THETAO , GA , GB , RO , EXTG 

Format statements for error messages. 

100 FORMAT(2X, 'Number of elements is more than can be handled.') 

150 FORMAT (2X, 'Number of boundary elements is more than 

1 can be handled. Ratio must be corrected.') 

200 FORMAT (2X, 'Number of boundary nodes is more than can be 
1 handled.') 

300 FORMAT ^2X, 'Number of nodes is more than can be handled.') 

400 FORMAT (2X, ' Band width is greater than can be handled.') 

500 FORMAT (2X, 'Number of interior nodes is more than can be 
1 handled,') 

600 FORMAT (2X, 'Size cf structure stiffness matrix is greater than 
1 can be handled.") 

700 FORMAT (//,2X, 'Execution aborted. Storage capacity exceeded.') 

800 FORMAT (2X f 'Ratio must be greater than or equal to 1 and less 

1 than or equal to 5.'# 

1 /,2x, 'Execution aborted.') 

Formats for output of information supplied by user. 


1000 FORMAT (// 2X, »>******»«**************** INPUT INFORMATION', 

1 ' »»**» ** * *»****«»*»* »*****' f /// , 

1 2X, 'MESH INFORMATION' ,/,2X, ' — — ' ,//, 


1 2X, 'SWEEP = ' ,F11. 6, 2X, 'NUMBER OF ANGLE DIVISIONS = ',15,//, 
1 2X, 'RANGE = ' ,F11. 6, 2X, 'NUMBER OF RADIAL DIVISIONS = ',15,//, 
1 2X,6X, 'SHAPE CONSTANT - Fll.fi, //, 

1 2X,16X,' INTERNAL BOUNDARY ELEMENTS',/, 

1 2X,6X, 'RATIO OF — =' ,15,/, 

1 2X,16X, 'EXTERNAL BOUNDARY ELEMENTS',//) 

1500 FORMAT (2X,6X, 'RADIAL DIVISIONS ARE GRADUATED',//) 

1600 FORMAT(2X,6X, 'RADIAL DIVISIONS ARE UNIFORM' ,/,29X, ' ' 

2000 FORMAT (2X, 'FREE FIELD WAVE INFORMATION' ,/,2X, — 

1 ' —',//, 


'• • *»• '*■*. • " t I ' t --• * *** *# 


1 2X FREQUENCY SQUARED * * ,F11.6,2X, 'FREQUENCY * »,F11.6,//, 
1 2X, 'ANGLE OF INCIDENCE FROM DOWNWARD VERTICAL * »,F11.6,/, 


1 2X, * (+ CLOCKWISE,- COUNTERCLOCKWISE) »,///, 

1 2X, 'MATERIAL PROPERTIES' ,/,2X, ' — — -',/A 


1 2X,' INTERNAL SHEAR MODULUS * »,F11.6,' Y + »,F11.6,//, 

1 2X,» INTERNAL MATERIAL DENSITY * ',FU.6,//, 

1 2X,' EXTERNAL SHEAR MODULUS -*,F11.6,//) 

C 

C Halt program execution, 

c 

STOP 

END 

C* ******** ******** ••*•»•••*•*******•*******•»****•*****•******************•**** 

c 

c«** ********************************** ***************************************** 
C SUBROUTINE MESH 
C — 

c 

C The purpose of the MESH subprogram is to generate the finite element 
C mesh required with a minimum of input by the user. It will generate a 
c mesh for a circular structure defined in polar- coordinates, 
c There are only four 

c required inputs. The virst is the angular sweep of the entire 
c structure in radians. The second is the maximum radius of the 
C structure. The two remaining inputs are simply the number of divisions 
C of these two structure dimensions that are desired. MESH will generate 
C two tables, The first is an incidence table for each element. The 
C second is a table of element coordinates and element type. These two 
C tables, GEOM and PROP are passed back to the main program as they are 
C used by a number of other subroutines in defining the stiffness 
C matrix, 
c 

SUBROUTINE MESH (PROP, GEOM) 

C 

e Declare in COMMON all variables required by more than one 
C subroutine. These variables are divided into four catagoriess 
C 1. SIZE This common block contains variables which define 

C array dimensions. 

C 2. PROB This common block contains variables which define 

C problem parameters. 

C 3. MAP This common block contains variables which define 

c the finite element mesh. 

c 4. dbug This common block contains variables which define 

C which debug switches are desired. 

c Arrays are not passed through COMMON but passed as arguments to the 
c subroutines, 
c 

COMMON /SIZE/ NUMANG , NUMRAD ,NtJMELM , INTNOD , NUMNOD , COD I. AG , 

1 NUMEQN , LENGTH .WIDTH , SPACE , BNDELM , BCOLS , ALENG , AWIDE, ATMPL 

INTEGER NUMANG , NUMRAD , NUMELM , INTNOD , NUMNOD , CODI AG , NUMEQN , 

1 LENGTH , WIDTH , SPACE , BNDELM , BCOLS , ALENG , AWIDE , ATMPL 

COMMON /PROB/ W2 , GA , GB , RO , EXTG , ALPHA , THETAO 
REAL W2 , GA , GB , RO, EXTG , ALPHA , THETAO 

COMMON /MAP/ SWEEP, RANGE, DRAD , DANG , BDANG1 , BDANG2 , RATIO , GRADNT 
INTEGER RATIO 

REAL SWEEP , RANGE , DRAD , DANG , BDANG1 , BDANG2 

LOGICAL GRADNT 

common /dbug/ dmain,dmesh,dassm,dktest,dforce, 

1 dhspac , drect , dtr i , dgamma , df ilam , damat , dbmat , 

1 dcmat , dnorm , ddiag , ddiag2 , df unc , dder iv , detest , datest 
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logical dmain , dmesh , das sm , dk test, df orce , dh spac , 

1 drect , dtri , dgamma , df ilam , damat , dbmat , dcmat , dnorm , 

1 ddiag ,ddiag2 ,df unc ,dderiv , detest ,datest 

Declare and dimension arrays required in this subroutine. 

INTEGER GEOM(NUMELM,4) 

REAL PROP (NUMELH #5) 

Declare variables required only within this subroutine. 

INTEGER I / J , ELM 
REAL RAD, THETA 

Calculate mesh sizes. 

DANG' = SWEEP/FLOAT (NUMANG) 

DRAD = RANGE/FLOAT (NUMRAD) 

BDANG1 - DANG*RATIO 

BDANG2 = SWEEP - (BNDELM-1) *BDANG1 

iterate over structure defining element coordinates and element-node 
incidences. 

Iterate over number of angle divisions. 

DO 20 I=0,NUMANG-1 

Define angle; o to sweep minus one angle mesh size. 
THETA = 1 1* DANG > 

Iterate over number of radius divisions, 

DO 10 J=0 , NUMRAD- 1 

Label element number. 

ELM=I * NUMRAD+ J+l 

Define element coordinates and types: 

1. if radial divisions are to be uniform or 

IF(GRADNT) GOTO 7 
PROP (ELM ,1) = ( RANGE- J* DRAD) -DRAD 
PROP ( ELM , 2 ) =RANGE- J * DR AD 
GOTO 8 

2. if radial divisions are to be graduated. 

7 PROP (ELM , 1) = ( (l-DANG/2 . 0) / (l+DANG/2.0) ) ** (0+1) 

1 "RANGE 

IK-7 . EQ . NUMRAD-1) PROP (ELM , 1 ) =0 . 0 
PROP ( ELM , 2 ) = ( ( l-DANG/2 . 0 ) / ( l+DANG/2 . 0 ) ) " • J 
1 "RANGE 

8 PROP ( ELM , 3 ) =THETA 
PROP ( ELM , 4 ) =THETA+DANG 
PROP (ELM, 5) =4.0 

IF (J.EQ. NUMRAD-1) PROP (ELM, 5) =3.0 
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J ’V 


C 

C Define element node incidences. 

C 

GEOM ( ELM , 1 ) ■ ( J+l ) * ( NUMANG+1 ) +1+1 

IF ( J . EQ . NUMRAD-1. ) GEOM (ELM , 1 ) “ (NUMANG+1 ) * NUMRAD+1 

GEOM (ELM,2)*J* (NUMANG+1 ) +1+1 

GEOM ( ELM i 3 ) «GEOM ( ELM , 2 ) +1 

GEOM (ELM, 4) “GEOM ( ELM , 1 ) +1 

IF ( J . EQ . NUMRAD-1) GEOM (ELK , 4) =0 

10 CONTINUE 

C 

C Check to see if mesh defines a complete circle. If so impose 
c constraint that coordinates and incidences along angle equal 0 and 
C 2Pi are the same. 

C 

IF(THETA.NE, (6. 28319-DANG) ) GOTO 20 
PROP (ELM, 3) “PROP (1,2) 

PROP (ELM, 4) “PROP (1,1) 

GEOM (ELM, 3) “GEOM (1,2) 

GEOM (ELM. 4) “GEOM (1,1) 

20 CONTINUE 

C++++++ ++++++++++++++++ +++++++++++++++++++++++++++++++++++++++++++++++++++ 

C If debug required output coordinate and incidence information. 

C 

if {.not. dmesh) return 
write (26, 110) 
do 30 i=l,numelm 

write (26, 100) (geom(i, j) , j=l,4) , (prop(i, j) , j=l,5) 

30 continue 

100 format(2x,4(i6) ,2x,5(fl0.3)) 

110 format (2x, 'nodel ,node2,node3 ,node4 ',2x,' ra',8x,' rb',8x, 

1 * qa',8x,' qb') 

C++++++ +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +++++++ 

c 

C Return to main program. 

C 

RETURN 

END 

C* ******************»******************, ********************* ******************* 

c* ***** * ************************************************************** ********* 

C SUBROUTINE ASSM 

C 

c 

c The purpose of the ASSM subprogram is to assemble the structure 
C stiffness matrix. It does this by looking at each element, obtaining 
c its element stiffness matrix and then placing it term by term into the 
C structure stiffness matrix. ASSM calls on the two subroutines RECT 
C and TRI to calculate element stiffness matrices, once the 
C structure stiffness matrix is assembled ASSM returns it and the 
c back to the main program for use in other subroutines. 

C Another vector passed form the main program to ASSM is the OUTPUT 
C vector. There is no information passed here other than the size of the 
c vector. The size of the output vector is dependent on the structure 
c size and therefore dimensioned in the main program for convenience. 

C ASSM will output the structure stiffness matrix if debug diagnostics 
c are requested. 

C 

SUBROUTINE ASSM (PROP , GEOM , K , OUTPUT ) 

C 
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C Declare in COMMON all variables required by more than ore 
C subroutine. These variables are divided into four catagoriess 
C 1. SIZE This common block contains variables which define 

C array dimensions. 

C 2. PROS This common block contains variables which define 

C problem parameters. 

C 3. MAP This coiii”-*'" block contains variables which define 

c the finite element mesh. 

C 4. dbug This common block contains variables which define 

c which debug switches are desired. 

C Arrays are not passed through COMMON but passed as arguments to 
C subroutines. 

C 


the 


c 


l 

l 


1 

l 

l 

l 


COMMON /SIZE/ NUMANG , NUMRAD , NUKELM , INTNOD , NUMNOD , CODIAG , 

NUM.EQN , LENGTH , WIDTH, SPACE , BNDELM , BCOLS , ALENG , AWIDE , ATMPL 
INTEGER NUMANG , NUMRAD , NUMELM ( I NTNOD , NUMNOD , CODIAG , NUMEQN , 

LENGTH , W I DTH , SPACE , BNDELM , BCOLS , ALENG , AW I DE , ATMPL 
COMMON /PROB/ W2 , GA , GB , RO , EKTG , ALPHA , THETAO 
REAL W2 , GA , GB , RO , EXTG , ALPHA , THETAO 

COMMON /MAP/ SWEEP , RANGE , DRAD , DANG , UDANG1 , BDANG2 , RATIO , GRADNT 
INTEGER RATIO 

REAL SWEEP , RANGE , DRAD , DANG , BDANG1 , BDANG2 

LOGICAL GRADNT 

common /dbug/ dmain,dmesh,dassm,dktest,dforce, 

dhspac , drect , dtri , dgamma , df ilam , damat , dbmat , 
dcmat,dnorm»ddiag,ddiag2 ,dfunc, dderiv, detest ,daiest 
logical dmain , dmesh , *assm, dktest , df orce , dhspac , 

drect , dtri , dgamma , df ilam , damat; , dbmat , demat , dnorm , 
ddiag , ddiag2 , dfunc , dderiv , detest , datest 


C Declare and dimension arrays required in this subroutine. 

C 

INTEGER GEOM (NUMELM, 4) 

REAL PROP (NUMELM, 5) ,K(NUMNOD, WIDTH) ,OUTPUT(NUMNOD) ,KELM(10) 

C 


c Declare variables required only in this subroutine. 

C 

INTEGER I , N , M , ELM , JTREF1 , JTREF2 , JTREF3 , JTREF4 , PULL , ROW , 

1 COL, PUSH 

C 

C Initialize structure stiffness matrix to zero. 

C 

DO 15 1=1, NUMNOD 

DO 16 J=l, WIDTH 

K(I,J)=0.0 

16 CONTINUE 

15 CONTINUE 
C 

C Iterate through each element calculating stiffness matrix and then 
C loading into structure stifffness matrix. 

C 

DO 50 ELM=1, NUMELM 

C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++4-+++++++ 
C If debug required output element number, 

c 

if (dassm) write (26, 120) elm 
120 format(2x, 'element=' ,i4) 

C+++++++++++++++++++++++++++++++++++++++++4-++++++++++++++++++++++++++ 
C Determine element type and calculate proper element 

C stiffness matrix. 


non non ooooo ooo 


I=ELM 

IF(PR0P(ELM,5).EQ,4.0) CALL RECT { PROP , KELM , I ) 
IF(PR0P(ELM,5) .EQ.3.0) CALL TRI (PROP, KELM, I) 

Determine structure-element node reference points. 

JTREF1 = GEOM(ELM,l)-l 
JTREF2 = GEOM(ELM,2)-2 
JTREF3 = GEOM(ELM,3)-3 
JTREF4 = GEOM (ELM *4) -4 

Iterate through each term of element stiffness matrix and 
position it in structure stiffness matrix. Only one term of 
each symmetric pair is considered. 


30 

40 

50 


DO 40 N=l, int (PROP (ELM, 5) ) 

DO 30 M=l,int (PROP (ELM #5) ) 

IF (M.GT.N) GOTO 40 
PULL = M +N* (N-l)/2 
IF(N.LE.1)R0W = JTREF1 
I F ( N . GT . 1 ) ROW = JTREF2 
I F ( N . GT . 2 ) ROW ~ JTREF3 
IF(N.GT.3)R0W = JTREF4 
IF(M.LE.1)C0L = JTREF1 
I F ( M . GT . 1 ) COL = JTREF2 
IF(M.GT.2)COL - JTREF3 
IF(M.GT.3)COL = JTREF4 

I=ROW- (l+CQBIAG) 

PUSH=COL-I 

K (ROW , PUSH ) =K (ROW , PUSH) +KELM ( PULL) 

IF ( ROW. EQ. COL) GOTO 30 
I=COL- (l+CODIAG) 

PUSH=ROW-I 

K (COL , PUSH ) =K (COL , PUSH ) +KELM ( PULL) 

CONTINUE 

CONTINUE 

CONTINUE 


Multiply each stiffness term by 2 to represent half circle problem. 

DO 51 1=1 , NUMNOD 

DO 52 J=l, WIDTH 

K(l,J)=2.0*K(I,J) 

52 CONTINUE 

51 CONTINUE 

J-++++++++T++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 

If debug required output structure stiffness matrix. 

if ( . not. dassm) return 
write(26,130) 

format (/,2x, ' **** structure stiffness matrix ****',/, 

2X , ' band storage mode ') 

write (26,100) (i,i=0, width) 
do 53 i=l,numnod 

write (26,200) i , (k (i , j) , j=l , width) 

continue 

format (/, i7, 100 (i8)) 
format (2x,i 5, 2x, 100 (f 6.2, 2x) ) 


130 


53 

100 

200 
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C+++++++’f++++++'f++++++++++++++++++++++++ , f+++++ - H‘++++++4'++++++++++++++++++ 

Return to main program. 

RETURN 
END 

C* * » * * * * M ****** *************** ********** ******************************** * 

C 

£•*•*** «««*«■* »n****«i»*«****it* *********************** ********************* ****** 

C SUBROUTINE RECT 

C '- 1 

C 

C The purpose of this subprogram is to calculate the element stiffness 
C matrix for a rectangular sector of a circle. The shape functions used 
C for the representation are linear in both radius and angle. The 
c details of the equations used here can be seen in the accompanying 
C thesis. The inputs required are the element coordinate information, 

C shear modulus and density of the structure, and the frequency of the 
C incoming wave. This data is supplied through the ASSM subprogram and 
c the element stiffness matrix KELM is returned, 
c 

SUBROUTINE RECT (PROP , KELM , ELM) 

C 

C Declare in COMMON all variables required by more than one 
C subroutine. These variables are divided into four catagories; 

C 1. SIZE This common block, contains variables which define 

C array dimensions. 

c 2. PROB This common block contains variables which define 

C problem parameters. 

C 3. MAP This common block contains variables which define 

C the finite element mesh. 

C 4. dbug This common block contains variables which define 

C which debug switches are desired. 

C Arrays are not passed through COMMON but passed as arguments to the 
c subroutines. 

C 

COMMON /SIZE/ NUMANG , NUMRAD , NUMELM , INTNOD , NUMNOD , CODIAG , 

1 NUMEQN , LENGTH , W IDTH , SPACE , BNDELM , BCOLS , ALENG , AW I DE , ATMPL 

I NTEGER NUMANG , NUMRAD , NUMELM , INTNOD , NUMNOD , CODIAG , NUMEQN , 

1 LENGTH , WIDTH, SPACE , BNDELM , BCOLS , ALENG , AW I DE , ATMPL 

COMMON /PROB/ W2 , GA , GB , RO , EXTG , ALPHA , THETA 0 
REAL W2 , GA , GB , RO , EXTG , ALPHA , THETAO 

COMMON /MAP/ SWEEP, RANGE, DRAD, DANG, BDANG1 , BDANG2 , RATIO, GRADNT 
INTEGER RATIO 

REAL SWEEP, RANGE, DRAD, DANG, BDANG1 , BDANG2 
LOGICAL GRADNT 

common /dbug/ dmain,dmesh,dassm,dktest,dforce, 

1 dhspac , drect , dtri ,dgamma , df ilam , damat , dbmat , 

1 dcmat ,dnorm,ddiag,ddiag2 ,dfunc, dderiv, detest , date st 

logical dmain,dmesh,dassm,dktest,dforce, dhspac, 

1 drect , dtri , dgamma , df ilam , damat , dbmat , dcmat , dnorm , 

1 ddiag , ddiag2 , df unc , dderiv , detest , datest 

C 

C Declare and dimension arrays required in this subroutine. 

C KELM, KSTIF, and KM ASS are element stiffness dependent and are therfore 
C dimensioned with numbers here as they would have to be changed if the 
C type of element were to be changed. They are independent of the 
C structure size. KELM must also be dimensioned in ASSM. 

C 

REAL PROP (NUMELM, 5) ,KELM(10) ,KSTIFF(10) ,KMASS(10) ,KSTIF2(10) 
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c 

c Declare variables required only by this subroutine. 

C 

INTEGER ELM 

REAL RA,RB,QA,QB, OIFFQ , D I FFR , RD I F$Q , DENOM 

C 

C Assign simplified variable names, 
c 

RA = PROP (ELM, 1) 

RB = PROP (ELM, 2) 

QA = PROP (ELM, 3) 

QB ~ PROP (ELM, 4) 

DIFFQ = QB-QA 
DIFFR = RB-RA 
RDIFSQ = RB* *2 - RA**2 
DENOM * 6 . 0*DIFFQ*DIFFR**2 
C$ IF (GA.EQ.O.O)GOTO 17 

C 

C Calculate element stiffness properties. Equations are derived in the 
c text. Stiffness is composed of a variable stiffness term, a constant 
c stiffness term, and a mass term. These 

C are calculated seperately and then combined later in this subroutine. 

C If debug is required each term is output. 

C 

c — ---------- 

KSTIF2 (l) =( (RB**3-RA*«3)/(3.0*DIFFQ*DIFFR) « (DIFFQ* *2*COS(QA) 

1 +2 . 0*DIFFQ*SIN (QA) -2.0* (COS (QA) -COS (QB) ) ) ) 

1 - ( (COS (QB) -COS (QA) ) / (3 . 0*DIFFQ) *DIFFR**2) 

G ++++++++++++++++++++++++++++++++ 

if (drect)write (26,141) 

141 format (2x, 'kstif 2 (1-10) ') 

if (drect)write(26,*)kstif2(l) 

C ++*f4*Hh+4'4-+*f*f 4* 4'+++++++++++ , f ++++ 

C 

KSTIF2 (2) =( (RB**3-RA**3)/(3 . 0*DIFFQ*DIFFR) * (- (DIFFQ**2) *COS (QA) 
1 -2 . 0*DIFFQ*SIN (QA) +2 . 0* (COS (QA) -COS (QB) ) ) ) 

1 - ( (COS (QB) -COS (QA) ) / (6. 0*DIFFQ) *DIFFR* *2) 

C ++++++++++++++++++++++++++++++++ 

if (drect) write (26, *)kstif2 (2) 
c ++++++++++++++++++++++++++++++++ 

c 

KSTIF2 (3) = ( (RB* *3-RA**3 ) / (3 . 0* DIFFQ* DIFFR) * (DIFFQ* *2*COS (QA) 

1 +2.0*DIFFQ*SIN(QA)-2.0*(COS(QA)-COS(QB))) ) 

1 - ( (COS (QB) -COS (QA) ) / (3 . 0* DIFFQ) *DIFFR**2) 

C ++++++++++++++++++++++++++++++++ 

if (drect) write (26,*)kstif2 (3) 

C -f+4 , +4**f ++*f+++*f 4-+4‘4‘4 , 4 , 4**f4'+«f+4*+4 , 4'+++ 

C — 

KSTIF2 (4) =((RB**3-RA**3)/(3.0*DIFFQ*DIFFR)* (DIFFQ* (SIN(QA) 

1 +SIN (QB) ) -2 . 0* (COS (QA)-COS (QB) ) ) ) 

1 + ( (COS (QB) -COS (QA) ) / (6. 0*DIFFQ) *DIFFR**2) 

if (drect)write (26,*)kstif2 (4) 

C +++++++++++++++++++++++++++t++++ 

C 

KSTIF2 (5) = ( (RB* *3-RA**3) / ( 3 . 0*DIFFQ* DIFFR) * (-DIFFQ* (SIN (QA) 

1 +SIN (QB) ) +2 . 0* (COS (QA) -COS (QB) ) ) ) 

1 + ( (COS (QB) -COS (QA) ) / (3 . 0* DIFFQ) *DIFFR**2) 

C ++++++++++++++++++++++++++++++++ 

if (drect)write(26,*)listif2(5) 
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++++++++++++++++++++++++++++++++ 


KSTIF2 ( 6 ) «((RB**3-RA«»3)/(3.0i»DIFFQ*DIFFR)» (-(DIFFQ**2)*C0S(QB) 
1 +2 . 0*DIFFQ*SIN (QB) -2 .0* (COS (QA) -COS (QB) ) ) ) 

1 -( (COS (QB) -COS (QA) )/ (3 .0*DIFFQ) »DIFFR«»2) 

+++++++++++++++++++++++++f++++++ 

if (drect) write (26,*) kstif 2 (6) 

+++++++ - +++-t+++++++++++++++++++++ 


KSTIF2 (7) = ( (RB**3-RA* *3)/ (3 .0*DIFFQ»DIFFR) * (-DIFFQ* (SIN (QA) 
1 +SIN (QB) ) +2 . 0* (COS (QA) -COS (QB) ) ) ) 

1 + ( (COS (QB) -COS (QA) ) / (3 . 0*DIFFQ) *DIFFR**2) 

++++++++++++++++++++++++-H+++++++ 
if (drect)write(26,*)ksfcif2(7) 
++++++++++++++++++++++++++++++++ 


KSTIF2 ( 8 ) =( (RB*»3-RA**3)/(3.0*DIFFQ*DIFFR)*(DIFFQ*(SIN(QA) 
1 +SIN (QB) ) -2 . 0* (COS (QA) -COS (QB) ) )) 

1 + ( (COS (QB) -COS (QA))/(6.0«DIFFQ)»DIFFR»*2) 

+++++++++++++++++++++++++-H-+++++ 

if (drect) write (26, *) kstif 2 (8) 
++++++++++++++++++++++++++++++++ 


KSTIF2 (9) = ( (RB* * 3-RA* *3) / (3 . 0*DlFFQ*DIFFR) * (DIFFQ**2*COS (QB) 
1 -2 . 0*DIFFQ*SIN (QB) +2.0* (COS (QA) -COS (QB) ) ) ) 

1 - ( (COS (QB) -COS (QA) )/ ( 6 . 0*DIFFQ) *DIFFR»*2) 

+++++++++++++ 4-++++++++++++++++++ 

if (drect)write (26, *)kstif2 (9) 

++++++++++++++++++++++++++++++++ 


KSTIF2 (10) = ( (RB* *3-RA* *3) / (3 . 0*DIFFQ*DIFFR) * (- (DIFFQ**2) »COS(QB) 
1 +2 .0*DIFFQ*SIN(QB) -2 .0« (COS (QA) -COS (QB) ) ) ) 

1 - ( (COS (QB) -COS (QA) ) / (3 .0*DIFFQ) «DIFFR**2) 

•f 4- *f + + 4* 4* + + 4* 4* 4* 4* 4* + + + + 4* + + + + + 4* 4* 4* 4* + + 4- + 
if (drect) write (26, *) kstif 2 (10) 

+++++ 4 -++++++++++++++++++++++++++ 


17 KSTIFF (1) = (DIFFQ**2*RDIFSQ + (-9 . 0*RB**2+12 . 0*RA*RB-3 . 0*RA* *2) + 
1 6.0*RB**2*ALOG(RB/RA))/DENOM 

4* 4* 4* 4- + 4* + 4* 4* 4* + + 4* + 4* + 4* + 4* 4 4* 4 4 4 4 4 4 4 4* 4 4 4 
if (drect) write (26,140) 

140 format (2x, 'kstif f (1-10) ') 

if (drect)Write(26,*)kstiff (l) 

+++++++++ H-+++++++++++++++++++++ 


KSTIFF (2) = (- (DIFFQ**2) *RDIFSQ + 3.0*RDIFSQ - 
1 6 . 0*RA*RB* ALOG (RB/RA) ) /DENOM 

++++++++++++++++++++++++++++++++ 
if (drect ) write (26 , *) kstif f (2 ) 
++++++++++++++++++++++++++++++++ 


KSTIFF (3) = (DIFFQ**2*RDIFSQ + (3.0*RB**2-12.0*RA*RB+9.0»RA«*2) + 
1 6.0*RA**2*ALOG(RB/RA))/DENOM 
C ++++++++++++++++++++++++++++++++ 

if (drect)write(26,*)kstiff (3) 

C ++++++++++++++++++++++++++++++++ 

C * 

KSTIFF (4) = (-0 .5*DIFFQ**2*RDIFSQ - 3.0*RDIFSQ + 

1 6.0*RA*RB*ALOG(RB/RA) ) /DENOM 
++++++++++++++++++++++++++++++++ 


C 


OOOQOQOO OQOQ OOO OOO OOO OOO OOO O O 


if (drect) write ( 26, *) kstif f (4) 
++++++++++++++++++++++++++++++++ 


KSTIFF(5) B (0.5«DIFFQ»«2"RDIFSQ + (-3 .0«RB**24-X2.0'»RA»RB~9.0*RA«»2) 
X -6.0*RA**2*AL0G(RB/RA) )/DENOM 
++++++++++++++++++++++++++++++++ 
if (drect) write (26,*) kstif f (5) 

+++++++++++4’^ +++++■++++ ■<"♦•+++++++++ 


KSTIFF (6)=(DIFFQ**2*RDIFSQ 4- (3 .0«RB«<*2-12 ,0*RA*RB+9.0*RA*»2) 4- 
1 6.0»RA»*2*ALOG(RB/RA))/DENOM 

4-++4-++++++++++4 , ++ - f++++++++++'f+++ 
if (drect ) write (26, *) kstif f (6) 

++++++++++++++++ , +4'+++’H"++++++4-++ 


KSTIFF (7) =(0.5*DIFFQ»«2*RDIFSQ 4- (9.0»RB*«2-X2.0'»RA»RB4-3.0«RA»«2) 
X - 6.0*RB*a 2«AL0G(RB/RA))/DEN0M 
++++++++++++++++++++++++++++++++ 
if (drect) write (26, *) kstif f (7) 

++++++++++++++++++++++++++++++++ 


KSTIFF(8)=(-0.5«DIFFQ»*2*RDIFSQ - 3.0»RDIFSfi + 
X 6.0*RA*RB*ALOG(RB/RA) )/DENOM 

++++++++++++++++++++++++++++++++ 
if (drect) write (26, ») kstif f (8) 

++++++++++++++++++++++++++++4"f++ 


KSTIFF(9)=(-(DIFFQ**2) *RDIFSQ + 3.0*RDIFSQ - 
X 6 . 0*RB*RA*ALQG (RB/RA) ) /DENOK 
++++++++++++++++++++++++++++++++ 
if (drect ) write (26, *) kstif f (9) 
++++++++++++++++++++++++++++++++ 


KSTIFF ( 10) = (DIFFQ**2*RDIFSQ + (-9.0«RB»*2+X2.0»RA*RB-3.0<»RA»*2) + 
X 6.0 «RB**2*ALOG(RB/RA) )/DENOM 
++++++++++++++++'{•+++++++++++++++ 
if (drect) write (26, *) kstif f (X0) 

++++++++++++++++++++++++++++++++ 


KMASS(X) =(DIFFQ) * (RB**2+2.0*RA*RB-3.0*RA**2)/36.0 
+ *f + + + + 4* + + + + *r *f + + + + + *f + + *f 4- + + + *f + + 
if (drect) write (26, X45) 

X45 format (2x, 1 kmass (X-XO) ' ) 

if (drect) write (26, * ) kmass (X) 
++++++++++++++++++++++++++++++++ 


KMASS(2)=(DIFFQ) *(RDIFSQ)/36.0 
++++++++++++++++++++++++++++++++ 
if (drect ) write (26, *>) kmass (2) 

++++++++++ +++4*++"l"++i : f*+++++++++4 , ++ 


KMASS(3)=(DIFF2)*(3.0*RB**2-2.0»RA*RB-RA«*2)/36.0 

■f +*f •f+++++‘f+4 , ’f*f *f +++++4 , ++ 

if (drect) write (26,*) kmass (3) 
++++++++++4-+++++++++++++++++++++ 


KMASS (4) = (DIFFQ) * (RDIFSQ)/72 . 0 
C 4- 4-4- 4" 4> 4- 4- 4- 4- 4-4 4- 4- 4- 4- 4- 4- 4- 4-4- 4- 4- 4- 4-4- 4- 4- 4- 4- 4-4- 4- 

if (drect ) write ( 26 , * ) kmass ( 4) 

4- 4^ 4- 4- 4- 4- 4- 4- 4- 4- 4- 4- 4- 4-4- 4-+4- 4-4- 4- 4^ 4- 4- 4- 4- 4- 4- 4- 4-4- 4- 


C 


on o ono no o on o on 


KMASS(5) = (DIFFQ)«(3.0*RB»*2-2.0*RA<'RB-RA**2)/72.0 
C ++++++++++++++++++++++++++++++++ 

if (drect) write (26, •) kmass (5) 

++++++++++++++++++ 4 +++++++++++++ 


KMASS (6) e (DIFFQ) * (3.0*RB**2-2.0*RA*RB-RA**2)/36.0 
++++++++++++++++++++++++-f +++++++ 
if (drect ) write (26,1*) kmas s (6) 

+++++++++++■♦• +++++++ 4 -+++++++++++ 4 - 


KMASS (7) = (DIFFQ) * (RB**2+2 . 0*RA*RB~3 .0*RA**2)/72.0 
++++++++++++++t++++t++++++++++++ 
if (drect) write (26, * ) kmass (7) 
++++++++++++++++++++++++++++++++ 


KMASS (8) -(DIFFQ) * (RDIFSQ)/72.0 
++++++++++++++++++++++++++++++++ 
if (drect) wr*.te (26, * ) kmass (8) 
++++++++++++++++++++++++++++++++ 


KMASS (9)«(DIFFQ) * (RDIFSQ) /36 . 0 
++++++++++++ 4 -+++++++++++++++++++ 
if (drect)wrlte(26,*)kmass(9) 
++++++++++++++++++++++++++++++++ 


KMASS (10) =(DIFFQ) * (RB**2+2 . 0*RA*RB-3 . 0*RA* *2) /36 . 0 
C +*f+*f i f++++4-+++++4*4*4-++++ , f4*4-f+4’+*f++ 

if (drect)write(26,*)kmass(lo) 
c ++++++++++++++++++++++++++++++++■ 

c 

C Combine stiffness and mass terms. 

C 

KELM ( 1 ) =RO * W2 * KMASS ( 1 ) -GB*KSTIFF (1) -GA*KSTIF2 (l) 

KELM (2) =RO*W2*KMASS (2) -GB*KSTIFF (2) -GA*KSTIF2 (2) 

KELM (3) =RO*W2*KMASS (3) -GB*KST IFF (3) -GA*KSTIF2 (3) 

KELM (4)=RO*W2*KMASS (4) -GB«KSTIFF (4) -GA*KSTIF2 (4) 

KELM ( 5) =R0*W2*KMASS (5) -GB*KSTIFF (5) -GA*KSTIF2 (5) 

KELM (6)=RO*W2*KMASS (6) -GB*KSTIFF (6) -GA*KSTIF2 (6) 

KELM (7) =RO*W2* KMASS (7) -GB*KSTIFF (7) -GA*KSTIF2 (7) 

KELM (8) =RO*W2*KMASS (8) -GB*KSTIFF (8) -GA*KSTIF2 (8) 

KELM (9) =RO*W2* KMASS (9) -GB*KSTIFF (9) -GA*KSTIF2 (9) 

KELM (10) =RO*W2*KMASS (10) -GB*KSTIFF ( 10) -GA*KSTIF2 (10) 

C 

C Return to, ASSM subroutine, 

c 

RETURN 

END 

C************************************************************ ft* ******* ********* 

c 

C*************************«**************K***********************ft*ft*********** 

C SUBROUTINE TRI 

C 

C 

C The purpose of this subprogram is to calculate the element stiffness 
c matrix for a triangular sector of a circle. The shape functions used 
C for the representation are linear in both radius and angle. The 
C details of the equations used here can be seen in the accompanying 
C thesis. The inputs required are the element coordinate information, 

C shear modulus and density of the structure, and the frequency of the 
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C incoming wave. This data is supplied through the ASSM subprogram and 
C the element stiffness matrix KELM is returned. 

C 

SUBROUTINE TRI (PROP, KELM, ELM) 

C 

C Declare in COMMON all variables required by more than one 
C subroutine. These variables are divided into four catagoriest 
C 1, SIZE This common block contains variables which define 

C array dimensions. 

c 2. PROB This common block contains variables which define 

C problem parameters . 

C 3. MAP This common block contains variables which define 

C the finite element mesh. 

C 4. dbug This common block contains variables which define 

C which debug switches are desired. 

C Arrays are not passed through COMMON but passed as arguments to the 
C subroutines. 

C 

COMMON 

1 

INTEGER 

1 

COMMON 
REAL 
COMMON 
INTEGER 
REAL 
LOGICAL 
common 

1 
1 

logical 

1 
1 
c 

C Declare and dimension arrays required in this subroutine. 

C KELM, KSTIF, and KMASS are element stiffness dependent and are therfore 
C dimensioned with numbers here as they would have to be changed if the 
C type of element were to be changed. They are independent of the 
C structure size. KELM must also be dimensioned in ASSM. 
c 

REAL PROP (NUMELM, 5) , KELM (10) ,KSTIFF(10) ,KMASS(10) ,KSTIF2(10) 
Declare variables required only by this subroutine. 

INTEGER ELM 

REAL RA,RB,QA,QB,DIFFQ 

Assign simplified variable names. 

RA = PROP (ELM, 1) 

RB = PROP (ELM, 2) 

QA = PROP ( ELM , 3 ) 

QB = PROP (ELM ,4) 

DIFFQ = QB-QA 
C$ IF(GA.EQ.O.O) GOTO 17 

C 

C Calculate element stiffness properties. Equations are derived in the 
C text. Stiffness is composed of a variable stiffness term, a constant 
C striffness term, and a mass term. These 


/SIZE'/ NUMANG , NUMRAD , NUMELM , INTNOD , NUMNOD , CODI AG , 

NUMEQN , LENGTH , WIDTH , SPACE , BND,ELM , BCOLS , ALENG , AW IDE , ATM PL 
NUMANG , NUMRAD , NUMELM , INTNOD , NUMNOD , CODI AG , NUMEQN , 

LENGTH , WIDTH , SPACE , BNDELM , BCOLS , ALENG , AWIDE , ATM PL 
/PROB/ W2 , GA , GB , RO , EXTG , ALPHA , THETAO 
W2 , GA , GB , RO , EXTG , ALPHA , THETAO 

/MAP/ SWEEP , RANGE , DRAD , DANG , BDANG1 , BDANG2 , RATIO , GRADNT 
RATIO 

SWEEP , RANGE , DRAD , DANG , BDANG1 , BDANG2 
GRADNT 

/dbug/ dmain,dmesh,dassm,dktest,dforce, 
dhspac,drect,dtri,dgamma,dfilam,damat,dbmat, 
dcmat , dnorm , ddiag , ddiag2 ,df unc ,dderiv , detest , datest 
dmain , dmesh , dassm , dktest , df orce , dhspac , 
drect , dtri , dgamma , df ilam ,damat , dbmat , dcmat , dnorm , 
ddiag, ddiag2,df unc, dderiv, detest , datest 
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are calculated seperately and then combined later in this subroutine. 
If debug is required each term i3 output. 


kUfIF2 (1)“ RB**2/ (3.0*DIFFQ) * {- (DIFFQ* »2) * (COS(QB)-COS(QA) ) ) 
++++++++++++++++++++++++++++++++■++ 
if (dtri)write(26,151) 

151 format (2x,')cstif 2 (1-6)*) 

if (dtri) write (26, »)kstif 2(1) 
++++++++++4 , ++++++++ , t , +++++t++++++++ 


KST1F2 (2) = RB* *2/ (3 . 0*DIFFQ) * (- (DIFFQ* *2) *COS (QA) +DIFFQ* 
1 (SIN (QB) -SIN (QA) ) ) 

++++++++++++++++++++++++++++++++++ 
if (dtri)write(26,«)kstif2(2) 
++++++++++++++++++++++++++++++++++ 


KSTIF2 (3)= RB**2/ (3 . Q*DIFFQ) » ( (DIFFQ* *2) *COS (QA) +■ * 
1 2 > 0*DIFFQ*SIN (QA) +1 , 0* (COS (QB) -COS (QA) ) ) 

+++-f++++++++++++++++++++++++++++++ 

if (dtri)write(26,*)kstif2(3) 
++++++++++++++++++++++++++++++++++ 


KSTIF2 (4)- RB* *2/ (3 . 0*DIFFQ) * ( (DIFFQ* *2) *COS (QB) -DIFFQ* 
1 (SIN (QB)-SIN (QA) ) ) 

+ + + + 4* *f + *f"f + + + *♦* •+” + + *f •f h ^ 4* + 4* + + *f + •+* + + 4* *f + 
if (dtri)write(26, *) Xstif2 (4) 

++++++++++++t+++++++++++++++++++++ 


KSTIF2 (5)- RB**2/(3.0*DIFFQ)«(-(DIFFQ**2.)* (SIN(QA)+SIN(QB) ) 
1 - (COS (QB) -COS (QA) ) ) 

++++++++++++++++++++++++++++++++++ 
if (dtri)write (26, « ) Xstif2 (5) 

+++4"+-f+++++++4*++4-4 < 4'4*‘f+'f*f+++++4' , f++4* 


KSTIF2 (6)= RB**2/(3.0*DIFFQ) * (DIFFQ**2*COS (QB) 

1 +2 . 0*DIFFQ*SIN (QB) + (COS (QB) -COS (QA) ) ) 

++++++++++++++++++++++++++++++++++ 
if (dtri) write (26, * ) Kstif 2 (6) 
++++++++++++++++++++++++++++++++++ 


17 KSTIFF (1) = 0.5*DIFFQ 

++++++++++++++++++++++++++++++++++ 
if (dtri) write (26, 150) 

150 format(2x, ’kstiff (1-6) 1 ) 

if (dtri) write (26, * ) kstif f (1) 
++++++++++++++++++++++++++++++++++ 


KSTIFF(2)= -0 . 25*DlFFQ 
++++++++-H-++++++++++++++++++++++-I-+ 
if (dtri) write (26, *) kstif f (2) 
++++++++++++++++++++++++++++++++++ 


KSTIFF (3)= (DIFFQ**2 + 3.0)/(6.0 * DIFFQ) 
++++++++++++++++++++++++++++++++++ 
if (dtri) write (26, *) kstif -f (3) 
++++++++++++++++++++++++++++++++++ 


KSTIFF(4)= -0.25*DIFFQ 
C ++++++++++++++++++++++++++++++++++ 


non qooo o on o on o no o oo o no o oo o no n no 


if (dtri) write (26,*) kstif f (4) 
+++++++++++++'f+4"f +++++++++++++++++ 


KSTIFF(5)* (DIFFQ* *2 - 6.0)/(12.0 • DIFFQ) 
•f++4‘4 > ++4 , +++++’* , ++++4"H"f ++++++++++++ 
if (dtri) write (26, « ) Kstif f (5) 

+■ ++++++++.+++++++++++++++++++++++++ 


KSTIFF(6)» (DIFFQ**2 + 3.0)/(6.0 * DIFFQ) 
++++++4"t4 - ++++++++++++++++'f ++++++++ 
if (dtri) write (26, * ) Kstif f (6) 

++++++ 4 •+++++++++++++++++++■♦• ++++++ + 


KMASS (1)« (DIFFQ) *RB**2/l2.0 
++++++++++++++++++++++++++++++++++ 
if (dtri)write(26,160) 

160 format (2x, 'kmass (1-6) ' ) 

if (dtri) write (26,*) kmass (l) 
++++++++++++++++++++++++++++++++++ 


KMASS (2) = (DIFFQ) *RB« *2/24.0 

++++++++++4'4 , +4'++++++<++++++++++++++ 
if (dtri) write (26, •) kmass (2) 

•f *f 4* 4* + 4* + *f + + 4* «f 4* 4* 4* + Hh 4- 4* + 4* 4* + 4* + 4* + 4* 4* 4* 4* 4* 4* 4* 


KMASS (3 )= (DIFFQ) *RB**2/l2.0 

+++++++++++++++++++++++++++■+++++++ 
if (dtri) write (26, * ) Kmass ( 3 ) 
++++++++++++++++++++++++++■++++++++ 


KMASS (4) = (DIFFQ) *RB* #2/24. 0 

+++++++4 - ++ i f+++++++++++++++++++++++ 
if (dtri) write (26, *) Kmass (4) 
++++++++++++++++++++++++++++++++++ 


KMASS (5 ) = (DIFFQ) *RB** 2/24.0 

++++++++++++++++++++++++++++T+++++ 

if (dtri) write (26, * ) kmass (5) 
++++++++++-(-+++++++++++++++++++++++ 


KMASS (6) = (DIFFQ) *RB*« 2/12.0 

+++++++++++++++++++++++'+++++++++++ 
if (dtri) write (26, * ) Kmass (6) 
+++++++++++♦++++++++++++++++++++++ 


Combine stiffness and mass terms. 

KELM (1) =R0*W2*KMASS (1) -GB*KSTIFF (l) -GA*KSTIF2 (l) 
KELM (2) =RO*W2*KMASS (2) -GB*KSTIFF (2) -GA*KSTIF2 (2) 
KELM ( 3 ) =RO*W2*KMASS ( 3 ) -GB*KSTIFF ( 3) -GA*KSTIF2 (3 ) 
KELM (4)=RQ*W2*KMASS (4)-GB»KSTIFF (4) -GA*KSTIF2 (4) 
KELM (5) =RO*W2«KMASS (5) -GB’*KSTIFF (5) -GA*KSTIF2 (5) 
KELM (6)=RO*W2*KMASS (6) -GB*KSTIFF (6) -GA*KSTIF2 (6) 

Return to ASSM subroutine. 

RETURN 

END 


o o o non 


C* »**»**«•* (m»l» • » »»**** Ik *.##* **#*»» ****** ** ** *»*»*»» M ■»*»»»»«»« *» 

SUBROUTINE KTEST (K, WO , WG , KOG , B , KOO , XL , OUTPUT , EXACT) 

C 

C Declare in COMMON all variables required by more than one 
C subroutine. These variables are divided into four catagoriess 
C 1. SIZE This common block contains variables which define 

c array dimensions. 

c 2. PROS This common block contains variables which define 

c problem parameters. 

C 3. MAP This common block contains variables which define 

C the finite element mesh. 

c 4. dbug This common block contains variables which define 

c which debug switches are desired. 

C Arrays are not passed through COMMON but passed as arguments to the 
C subroutines. 

C 

COMMON /SIZE/ NUMANG , NUMRAD , NUMELM , INTNOD , NUMNOD , CODI AG , 

1 NUMEQN , LENGTH /WIDTH , SPACE , BNDELM , BCOLS , ALENG , AWIDE , ATMPL 

INTEGER NUMANG , NUMRAD , NUMELM , INTNOD , NUMNOD , CODI AG , NUMEQN , 

1 LENGTH .WIDTH .SPACE , BNDELM , BCOLS , ALENG , AW IDE , ATMPL 

COMMON /PROB/ W2 . GA , GB . RO , EXTG . ALPHA > THETAO 
REAL W2 , GA . GB , RO , EXTG , ALPHA . THETAO 

COMMON /MAP/ SWEEP, RANGE, DRAD, DANG, BDANG1.BDANG2, RATIO, GRADNT 
INTEGER RATIO 

REAL SWEEP , RANGE , DRAD , DANG , BDANG1 , BDANG2 
LOGICAL GRADNT 

common /dbug/ dmain,dmesh,dassm,dktest,dforce, 

1 d h spac , d rfect , dt ri , dgamma , d f ilam , damat , dbmat , 

1 dcmat,dnorm,ddiag,ddiag2,dfunc,dderiv, detest, datest 

logical dmain,dmesh,dassm,dktest,dforce,dhspac, 

1 drect , dtri , dgamma , df ilam , damat , dbnk'it , demat , dnorm , 

1 ddiag,ddiag2,dfunc,dderiv, detest, datest 

Declare and dimension arrays required by this subroutine. 

REAL K (NUMNOD, WIDTH) ,WO(NUMEQN) .WG(INTNOD) ,KOG(NUMEQN, INTNOD) 

1 ,B (NUMEQN, INTNOD) , KOO (NUMEQN, WIDTH) , XL (NUMEQN , SPACE) , 

2 OUTPUT (NUMNOD) , EXACT (NUMEQN) 

Declare variables required by this subroutine. 

, INTEGER N.M.I, J , L, PULL, P , ROW, COL, IER 
C 

C Read in the values of displacements at the boundry nodes. These 
C values are computed from the free field motion at the boundary. 

C 

READ (25,0 (WG (I) ,1=1, INTNOD) 

C 

C Echo imposed boundary displacements. 

C 

WRITE(26, 725) 

DO 5 J=l, INTNOD 

WRITE(26,750) J,WG(J) 

5 CONTINUE 
C 

c Partition K( omega, gamma) from structure stiffness matrix. KOG is 
C pulled term by term from K. KOG is stored in full storage mode, 
c 

N=0 

DO 20 I=INTN0D+1, NUMNOD 
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N-N+l 

DO 10 J®1» INTNOD 

KOG(N,J)«0.0 
B(N, J)*Q.0 
M*I-(l+CODIAG) 

PULL«J»M 

IF(PULL.LE.O) GOTO XO 
IF (PULL. GT. WIDTH) GOTO 20 
KOG(N,J)*K(l,PULL) 

B(N, J)=K(I»PULL) »-X.O 

XO CONTINUE 

20 CONTINUE 
C 

c Partition K (omega, omega) from structure stiffness matrix. KOO is 
c pulled term by term from K. KOO is stored in band storage mode rather 
c than symmetric storage mode to aXXow the use of an IMSL solving 
c routine, 
c 

DO 40 I*lNTNOD+X,NUMNOD 

N®!- (INTNOD+I+CODIAG) 

P=N 

IF(N.LT.0)N“0 

H=I+CODIAG 

IF (M . GT . NUMNOD) M=NUMNOD 
DO 30 J=INTN0D+1+N,M 
L=I- (l+CODIAG) 

PULL=J-L 

ROW" I- INTNOD 

COL=J-INTNOD-P 

KOO {ROW , COL) =K ( I , PULL) 

30 CONTINUE 

40 CONTINUE 

C++++++++ - f+++++++++++++++++++-+++++++++++++++++++++++++++4+++++'f++++' , f +++•♦•++++++ 
If debug required output K (omega, gamma) and K (omega, omega) . 

if (.not.dktest) goto 500 
write(26,300) 

300 format(//,2x, ' *•** k (omega, gamma) »»**',/) 
do 50 i=l,numeqn 
do 60 j=X, INTNOD 
output (j)=kog(i, j) 

60 continue 

write (26*100) (output (n) ,n=X, INTNOD) 

50 continue 

write (26, 3X0) 

310 format (//,2x, ’ **** k( omega, omega) »»»•',/) 
do 70 i=l,numeqn 
do 80 j=l, width 
output ( j ) =koo (i , j ) 

80 continue 

write (26,200) (output (n) ,n=l, width) 

70 continue 

•V+++++++4.++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 

Solve equation KOO x X = B = --K0G. X, the solution is written over 
top of B. This solving is done with IMSL routine LEQT1B. 

500 CALL LEQT1B ( KOO , NUMEQN , CODI AG , COD I AG , NUMEQN , B , I NTNOD , NUMEQN , 0 , 

X XL, IER) 

IF(lER.EQ.129)write(5,400) 


52 


400 format (2x, * error ' ) 

C+++++++++++++++++++++++++++++++++++++++++++++++++++++ - f+’f+++4-f +++++++++++++++ 

C If debug required output irtSL altered KOO and solution matrix. 

C 

if (. not. dKtest) goto 600 
write (26, 320) 

320 format(//,2x, ' **** modified k( omega, omega) ***««,/) 
do 75 i=l,numegn 
do 85 j=l, width 
output (;J)=koo(i,;j) 

85 continue 

write (26,200) (output(n) ,n=l, width) 

75 continue 

write (26,330) 

330 format (//, 2x, ' **** solution matrix »***') 

do 90 i=l,numeqn 
do 95 j=l, INTNOTj 
output(j)*B(i, •}) 

95 continue 

writs (26,100) (output (n) ,n=l,INTN0D) 

90 continue 
100 format (2x,lllf6.2) 

200 format (2x,lllf6.2) 

C++++++++++++++-+++++++++++++++++++++++-H-+++++4-+++++++++++++++++++++++++++++++++ 

C 

C Calculate interior node displacements which equals the boundry node 
C displacements times the- solution matrix, 
c 

600 DO 110 I=1,NUMEQN 
W0(l)=0.0 
DO 120 J=l, INTNOD 

WO ( I ) =B (I , J) *WG (J) +WO (I) 

120 CONTINUE 

110 CONTINUE 
C 

C Iterate over nodal points solving for the real portion of the 
C displacements from the incident wave equation. This provides a 
C partial check on the correctness of the stiffness matrix. This check 
C is only valid for a uniform mass and shear modulus equal to those of 
C the exterior region. Thus the solved for solution is the free field 
C motion. 

C 

DO 25 I=0,NUMANG 

THETA = I * DANG 
DO 15 J=1,NUMRAD 

RAD = ( (l-DANG/2 . 0) / (l+DANG/2 . 0) ) **J* RANGE 
IF ( . NOT . GRADNT) RAD=RANGE- J * DRAD 
I F ( J . EQ . NUMRAD ) RAD=0 . 0 
N= ( J-l) * INTNOD+I+1 

EXACT (N)=0.5* (COS (SQRT (W2) *RAD*SIN (THETAD+THETAO) ) + 

1 COS (SQRT (W2 ) *RAD*SIN (THETA-THETAO) ) ) 

15 CONTINUE 

25 CONTINUE 
C 

C Output finite element solution and exact solution side by side for 
C comparison. 

C 

I F (. NOT . DKTEST ) RETURN 
WRITE (26, 700) 

DO 35 I=1,NUMEQN 


GOO OOO OOOOOOOGOOOQGO OOO OOO 


J-I+INTNOD 

WRITE (26# 750) J,W0(I) ,EXACT(I) 

35 CONTINUE 

700 FORMAT (// , 15X, * SOLVED * # 15X , ' EXACT ' ) 

725 FORMAT {//,10X,’ IMPOSED BOUNDARY DISPLACEMENTS PER',/# 
1 10X,* FREE FIELD MOTION ») 

750 FORMAT(2X,i5,3x,G15.6,6X,G15.0) 

Return to main program. 

RETURN 

END 


• «««*»**»i»»(l«K««»*«r*****»**»*««»«>«*«*f*»*»r****«********»**»**l**M*** 
SUBROUTINE FORCE(FG,FPHl) 

Declare in COMMON all variables required by more than one 
subroutine. These variables are divided into four catagories: 

1. SIZE This common block contains variables which define 

array dimensions. 

2. PROB This common block contains variables which define 

problem parameters. 

3. MAP This common block contains variables which define 

the finite element rc?sh. 

4. dbug This common block contains variables which define 

which debug switches are desired. 

Arrays are not passed through COMMON but passed as arguments to the 


subroutines. 


* 


COMMON 

/SIZE/ NUMANG, NUMRAD , NUMELM , INTNOD , NUMNOD , CODIAG , 

t. 

1 

INTEGER 

NUMEQN , LENGTH , WIDTH , SPACE, BNDELM , BCOLS , ALE?.{<3, AWIDE , ATMPL 
NUMANG , NUMRAD , NUMELM , INTNOD , NUMNOD #CODIAG , NUMEQN , 

1 

1 

COMMON 

LENGTH , WIDTH, SPACE , BNDELM , BCOLS , ALENG , AWI DE , ATMPL 
/PROB/ W2 , GA , GB , RO , EXTG , ALPHA , THETAO 

{i 

| 


REAL 

W2 , GA , GB , RO , EXTG , ALPHA , THETAO 

• l 


COMMON 

/MAP/ SWEEP , RANGE , DRAD , DANG , BDANG1 , BDANG2 , RATIO , GRADNT 



INTEGER 

RATIO 

1 


REAL 

SWEEP , RANGE , DRAD , DANG , BDANG1 , BDANG2 

1 

1 


LOGICAL 

GRADNT 


common 

/dbug/ dmain,dmesh,daf4!m,dktest,dforce, 

| 

1 


dhspac ,drect,dtri, dgamma , df ilam ,damat , dbmat , 


1 


dcmat , dnorm , ddiag , ddiag2 , df unc , dderiv , detest , date st 



logical 

dmain , dme sh , dassm , dkte st , df orce , dhspac , 


1 


drect , dt r i , dgamma , df ilam , damat , dbmat , dcmat , dnorm , 


1 


ddiag , ddiag2 , df Unc , dderiv , detest , datest 



Declare and dimension arrays used in this subroutine. 

COMPLEX FG(INTNOD) ,FPHI (BNDELMJ 
REAL XI (50) , WGT (50) 

Declare variable required only by this subroutine. 

INTEGER NUMPTS # 1 #ELM# J 
REAL KAPPA , KAPPAR , QA , QB , THETA 

COMPLEX FGCNST , FPCNST , FG1 , FG2 , FPHI 1 , NFG1 , NFG2 , NFPHI1 
C 

C Define required constants. 

C 



54 


O Q 


KAPPA = SQRT(W2*RANGE/EXTG) 

KAPPAR = KAPPA • RANGE 
FGCNST = CMPLX( 0.0, -KAPPA/4.0) 

FPCNST = CMPLX(0. 25,0.0) 

C ■++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
C If debug required output computed constants. 

C 

if (dforce) write (26, 1900) fgcnst,fpcnst 

C ++++++++++++++++++++++++++++++-f++++++++++++++++++++++++++++++++4-++++++++++ 
C 

C Read in Gauss integration points for numerical Integration of force 
C expressions, 
c 

READ(24, *) NUMPTS 
DO 30 I=l,NUMPTS/2 

READ (24, *) XI (I) ,WGT(I) 

30 CONTINUE 

c 

c 

c Compute forces at boundary nodes. 

C 

C 

c initialize forces. 

C 

DO 10 1=1, INTNOD 

FG(I)=(0. 0,0.0) 

10 CONTINUE 

C ++++++++++++++++-H-++++++ 

if (dforce) write (26, 1300) 

C ++++++++++++++++++++++++ 

C 

C Integrate over each element for both terms of shape function, 
c Contributions to common nodes are summed, 
c 

DO 20 ELM=1,NUMANG 

C ++++ 4*+++*f++*f +4'*f++4**f++Hh++'f , ++*f -f ++*f 4*+*f +*f*f*f 

j=elm 

if (dforce. and. elm. le.l) write (26, 300) elm 
c ■ +++++++++++++++++++++++++++++++++++++++ 

c 

C Define integration limits. 

C 

QA = (ELM-1) "DANG 
QB = ELM* DANG 

C ' +++++++++++++++++++++++++++++++++++++++++ 

if (dforce. and. elm. le.l) write ( 26, 400) qa,qb 
c +++++++++++++++++++++++++++++++++++++++++ 

C 

C Iterate through each Gauss integration point. 

C 

DO 50 1=1 , NUMPTS/2 
C 

C Transform -1 to 1 positive integration points to qa 

C to qb points. 

C 

THETA = ( (QB-QA) *XI ( I ) + (QB+QA) ) /2 . 0 
C +++++++++++++++++++++++++++++++++++++++++ 

if (dforce . and , elm . le . 1 ) write ( 26 , 500 ) theta 
+++++++++++++++++++++++++++++++++++++++++ 


ana 


Evaluate force function at positive integration 
point. 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


1 

c 

50 

20 


CALL GAMMA (£A , QB , THETA , KAPPAR , FG1 , FG2 , j ) 
+++++++++++++++++++++++++++++++++++++++++ 
if (af orce . ana . elm . le . 1) write (26, 700) fgl , f g2 

++++++++++i'’t++++++++++++++++4 , +++++ , f ++++++ 

Transform -1 to 1 negative integration points to ga 
to qb points, 

THETA = ((QA-CB)*XI(l)+(QB+QA))/2.0 
+++++++++++++++++++++++++++++++++++++++++ 
if (af orce .ana .elm . le .1) write (26,500) theta 

+++++++++++++++4++++++++++++++++++ +++++++ 

Evaluate force function at negative integration 
point. 

CALL GAMMA (QA , QB , THETA , KAPPAR , HPG1 , NPG2 , j ) 
+++++4-++++++++++++++-H-++++++++++++++-H-+++++++ 
if (af orce. ana. elm. Ie.l) write (26, 800) nfgl,nfg2 
+++++++++++++++++++++++++++++++++++++++++++++ 

Multiply evaluation by integration weight ana 
function constant ana sum. 

FG ( ELM ) =FG ( ELM ) +FGCNST* WGT (l)/2.0* ( FG1+NFG1 ) 

FG ( ELM+ 1 ) =FG ( ELM+1 ) +FGCNST * WGT (I)/2.0« (FG2+HFG2 ) 
+++++++++++++++++++++++++-+++++++++++++++++++++++ 
< *(af orce. ana. elm. ie.l) j=elm+l 
if (af orce. ana. elm. le. 1) write (26, 600 ) elm, fg (elm) , 
j,fg(j) 

++++++•*•+■++++++++++++++++++++++++++++++++++++++++ 

CONTINUE 

CONTINUE 


C — 

C 

C Compute forces over elements originating from the half space. 

c 

c 

C Initialize forces, 

c 

DO 80 I=1,BNDELM 

FPHI(I) = CMPLX(0. 0,0.0) 

80 CONTINUE 

C ++++++++++++++++++++-M-++ 

if (aforce)write (26,1400) 
c ++++++++++++++++++++++++ 

c 

C Integrate over each element, 

c 

DO 60 ELM=1,BNDELM 

C +++++++++++++++++++++++++++++++++++++++ 

j=elm 

if (af orce. ana. elm. ie.l) write (26, 300) elm 

c +++++++++++++++++++++++++++++++++++++++ 

c 

C Define integration limits. 

C 
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c 

c 

c 

c 

r, 

c 

c 

c 


c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


c 

c 

c 

c 


QA » (ELM-1) "BDANG1 
QB = ELM»BDANG1 

IF ( ELM . EQ . BNDELM ) QA = (ELM-1) *BDANG2 
I F ( ELM . EQ . BNDELM ) QB «= ELM*BDANG2 
+++++++++++++++++++++++++++4-+++++++++++-f+ 
if (df orce. and. elm. le.l) write (26, 400) qa,qb 

+++++++++++++++++++t+++++++++++++++++++-(-+ 

Iterate through each Gauss integration point. 

DO 70 1-1 ,NUMPTS/2 

Transform -1 to l positive Integration points to qa 
to qb points. 

THETA = ( (QB-QA) *XI (I) + (QB+QA) ) /2 .0 
++++++++++++++++++++++++-(-++++++++-t-+++++++ 

if (dforce. and. elm. le.l) write (26, 500) theta 
++++++++++++++++++++++++++++++■+++++++++++ 

Evaluate force function at positive integration 
point. 

CALL PH I LAM ( THETA , KAPPAR , FPH I 1 , j ) 
+++++++++++++++++++++++++++++++++++++++•++ 
if (dforce. and. elm. le.l) write (26, 1000) fphil 

+++++++++++++++++++++++++-H-++++++++++++++ 

Transform -1 to 1 negative integration points to qa 
to qb points. 

THETA = ((QA-QB) *XI ( I) + (QB+QA) ) /2 . 0 
+++++++++++++++++++++++++++++++++++++++++ 
if (dforce. and. elm. le.l) write (26, 500) theta 
++++++++++++++-K+++++++++++++++++++++++++ 

Evaluate force function at negative integration 
point. 

CALL PH I LAM ( THETA, KAPPAR , NFPH I 1 , j ) 
++++++++++++++++++++++++++++++++++++++++■+++ 
if (df orce. and. elm. le.l) write (26, 1100) nfphil 
++++++++++++++++++++♦++++++++++++++++++++++ 

Multiply evaluation by integration weight and 
function constant and sum. 


FPHI (ELM) =FPHI (ELK) +FPCNST* (QB-QA) * WGT ( S ) 

/2 . 0* (FPHI1+NFPHI1) 

++++++++++++++++++++++++++++ 4 -+++++++++++++++++++ 
if (df orce . and . elm . le . 1 ) write ( 26 , 1200 ) elm , 
fphi(elm) 

++T'+++++++++++++++++++++++++++++++++++++++++++++ 


70 CONTINUE 

60 CONTINUE 


C 

c 

C Multiply forces by 2 to represent half circle. 
C 

DO 65 R0W=2 , INTNOD-1 


FGfROW) =2.0*FG(ROW) 

65 CONTINUE 

DO 66 ROW=l,BNDELM 

FPHI (ROW) - 2 . 0*FPHI (ROW) 

66 CONTINUE 

C++++++++++++-K+++++++++++++++++++++++++++++++++++4-+++++++++++++++++++++++ 
C If debug required output FG, FPHI 
C 

if ( . not . dforce) RETURN 
write (26, >13,0) numpts 
write (26, 100) 
do 40 i*l,intnod 

write (26, 200) i,fg(i) 

40 continue 

write(26,150) 
do 45 i=l,bndelm 

write (26, 200) i,f phi (i) 

45 continue 
return 
C 

C Debug output formats, 

c 

100 format (//,2x, * node* ,14x, *fg' ,/,13x, ’real' ,12x, 'imaginary* ,/) 

110 format (/,2x, 'number of gauss points = *,i5) 

150 format (//,2x, ' node* ,14x, ' fphi’ ,/,13x, ' real' ,12x, 

1 'imaginary',/) 

200 format(2x,i5,2x,2gl5.7) 

300 format (2x, ' element =',i5) 

400 format(2x, 'qa ='gl5.7,/,2x, 'qb=* ,gl5.7) 

500 format(2x, 'theta ='gl5.7) 

600 format (2x, ' fg( » ,il, ' ) =' ,2gl5 .7,/,2x, *fg( ' ,il, ' ) =',2gl5.7) 

700 format (2x, 'fgl =' ,2gl5.7,/,2x, ’fg2 =',2gl5.7) 

800 format (2x, 'nfgi =' ,2gi5.7,/,2x, ’nfg2 =’,2gl5.7) 

1000 format (2x, 'fphil =',2gl5.7) 

1100 format (2x, 'nfphil =',2gl5.7) 

1200 format (2x, 'fph,i( ! ,11* ' ) =',2gl5.7) 

1300 format (2x , ' ******** fg check *****«*«*') 

1400 format (///,2x, ' *»***««* fphi check «•*******') 

1900 format (2x, 'fgcnst =' ,2gl5.7,/, 

1 2x,'fpcnst = ',2gl5.7) 

END 

C» * «»«»« ***** t************************ ******** **************** ********** 

C 

C*t**NkM******MA***A***** ft **M** * * * * ********************************** 

S UBROUTINE GAMMA ( QA , QB , THETA , KAPPAR , FG1 , FG2 , e lm ) 

C 

C Declare in COMMON all variables required by more than one 
C subroutine. These variables are divided into four catagories: 

C 1. SIZE This common block contains variables which define 

C array dimensions. 

C 2. PROB This common block contains variables which define 

C problem parameters. 

C 3. MAP This common block contains variables which define 

C the finite element mesh. 

C 4. dbug This common block contains variables which define 

C which debug switches are desired. 

C Arrays are not passed through COMMON but passed as arguments to the 
C subroutines. 

C 

COMMON /SIZE/ NUMANG, NUMRAD , NUMELM , INTNOD , NUMNOD , COD I AG, 
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1 NUMEQN , LENGTH , WIDTH , SPACE , BNDELM , BCOLS , ALENG , AWI DE , ATMPL 

I NTEGER NUMANG , NUMRAD , NUMELH , I KTNOD , NUMNOD , CODI AG , NUMEQN , 

1 LENGTH , WIDTH, SPACE , BNDELM , BCOLS , ALENG , AW I DE , ATMPL 

COMMON /PROB/ W2 , GA , GB , RO , EXTG , ALPHA , THETAO 
REAL W2 , GA , GB , RO , EXTG , ALPHA , THETAO 

COMMON /MAP/ SWEEP, RANGE, DRAD, DANG, BDANG1,BDANG2, RATIO, GRADNT 
INTEGER RATIO 

REAL SWEEP , RANGE , DRAD , DANG , BDANG1 , BDANG2 
LOGICAL GRADNT 

common /dbvg/ dmain,dmesh,dassm,d)stest,dforce, 

1 dhspac , drect , dtri , dgamma , df ilam , damat , dbmat , 

1 dcmat , dnorm , dd iag , ddiag2 , d f unc , dderiv , detest , dates t 

logical dmain ,dmesh ,dassm ,dktest ,df orce , dhspac , 

1 drect , d tri , dgamma , df ilam , damat , dbmat , dcmat , dnorm , 

1 ddiag , ddiagi'i , df unc , dderiv , detest , datest 

Declare variables used only in this subroutine. 

INTEGER elm 

REAL QA,QB, THETA, KAPPAR,C,D 

DOUBLE PRECISION SINDIF,SINSUM,CSSUM,SSDIF,SSSUM,A,B 
COMPLEX FG1 , FG2 , FGAMAA , NGAMA1 , NGAMA2 
C 

c Evaluate force function for the two terms of the shape function, 
c 

SINDIF « S' I N ( THETA'-T K ETAO ) 

S INSUM « SIN (THETA^ THETAO) 

CSDIF * COS (KAPPA’i^- SINDIF) 

CSSUM = COS (KAPPAR*S INSUM) 

SSDIF = SIN (KAPPAR«SINDIF) 

SSSUM = SIN (KAPPAR»SINSUM) 

A = SINSUH*CSSUM-SINDIF*CSDIF 
B = SINSUM*SSSUM+SINDIF*SSDIF 
C = A 
0 = B 

FGAMAA = CMPLX (C ,D) 

NGAMA1 = CMPLX ((QB-THETA) ,0.0) 

NGAMA2 = CMPLX ( (THETA-QA) ,0.0) 

FG1 = NGAMA1« FGAMAA 
FG2 = NGAMA2* FGAMAA 

C +++++++++++++++++++++++++++++++++++++++++++"!’4'+++'H-+'4-++++ 

C If debug required write results of computations. 

C 

if (dgamma. and. elm. le.l) write (26, 100) sindif , sinsum,csdif , 

1 cssum, ssdif , sssum, f gamaa , ngamal , ngama2 , f gl , f g2 
C +++++++++++++++++++++++++++■,'■++++++++++++++++++++++++++++ 

C 

C Format statements for debug output. 

C 

100 format (//,2x, 'sindif =',gl5.7,/, 

2x,'sinsum =',gl5.7,/, • 

2x, ' csdif =' ,gl5.7,/, 

2x, 'cssum =’,gl5.7,/, 

2x, 'ssdif =' ,gl5.7,/, 

2x, 'sssum =',gl5.7,/, 

2x,'fgamaa =',2gl5.7,/, 

2x, 'ngamal =• ,2gl5.7,/, 

2x,'ngama2 =',2gl5.7,/, 

2x,'fgl =' ,2gl5.7,/ , 

2x,'fg2 =' ,2gl5.7 ,//) 
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Return to FORCE subroutine, 

RETURN 

END 

**** *A*»**l>**«»******ft*»»** 


SUBROUTINE PHILAM (THETA , KAPPAR , FI , elm) 


Declare In COMMON all variables required by more than one 
subroutine. These variables are divided into four catagories: 

1. SIZE This common blocK contains variables which define 

array dimensions. 

2. PROB This common block contains variables which define 

problem parameters. 

3. MAP This common block contains variables which define 

the finite element mesh. 

4. dbug This common block contains variables which define 

which debug switches are desired. 

Arrays are not passed through COMMON but passed as arguments to the 
subroutines, 


1 

1 


1 

1 

1 

1 


COMMON /SIZE/ NUMANG, NUMRAD, NUMELM, INTNOD.NUMNOD.CODIAO, 

NUMEQN , LENGTH .WIDTH, SPACE , BNDELM , BCOLS , ALENG , AWIDE , ATMPL 
INTEGER NUMANG , NUMRAD , NUMELM , I NTNOD , NUMNOD , CODI AG , NUMEQN , ' 

LENGTH .WIDTH , SPACE , BNDELM , BCOLS , ALENG , AWIDE , ATMPL 
COMMON /PROB/ W2 .GA.GB.RO, EXTG, ALPHA, THETAO 
REAL W2 , GA , GB , RO , EXTG , ALPHA , THETAO 

COMMON /MAP/ SWEEP, RANGE, DRAD, DANG, BDANG1.BDANG2, RATIO, GRADNT 
INTEGER RATIO 

REAL SWEEP , RANGE , DRAD , DANG , BDANG1 , BDANG2 

LOGICAL GRADNT 

common /dbug/ dmain,dmesh,dassm,dktest,dforce, 

dhspac ,drect,dtri, dgamma , df ilam , damat , dbmat , 
dcmat , dnorm , ddiag , ddiag2 , df unc , dderiv , detest , datest 
logical dmain , dme sh , dassm , dkte st , df orce , dhspac , 

drect , dt r i , dgamma , df ilam , damat , dbmat , dcmat , dnorm , 
ddiag , ddiag2 , df unc , dderiv , detest , datest 


Declare variables used only in this subroutine. 


INTEGER elm 

REAL THETA, KAPPAR, C,D 

DOUBLE PRECISION SINDIF.SINSUM, CSSUM, S3PIF, SSSUM, A, B 
COMPLEX FI 

Evaluate force function. 

SINDIF = SIN (THETA-THETAO) 

SINSUM = SIN (THETA+THETAO) 

CSDIF = COS (KAPPAR*SINDIF) 

CSSUM = COS (KAPPAR*S INSUM) 

SSDIF = SIN (KAPPAR*SINDIF) 

SSSUM = SIN (KAPPAR»SINSUM) 

A = CSSUM + CSDIF 
B = SSSUM - SSDIF 
C = A 
D = B 

FI - CMPLX(C.D) 
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+++++++++++++++++ - f+ , f+++++++++++4 r ++'f+4'++++++++++4'++++4 , +++ 

If debug required output computed values. 

if (df ilam. and. (elm. le.l. or. elm. eg. 3000) ) write (26,100) sindif , 
1 sinsum,csdif ,cssum,ssdif ,sssum,fl 

4++4-+++++++-f-h+++++++++++++++++++++++++++++++++4"f++++++++ 

Format statements for debug output. 

100 format (//, 2x, * sindif =',gl5.7,/, 

2x,'sinsum =’,915.7,/, 

2x,'csdif =’,gl5.7,/, 

2x , 'cssum =' ,gi5.7,/, 

2x, ’ ssdif = ! ,gl5.7,/, 

2x, ' sssum =’ ,gl5.7,/, 

2 x, ' £1 =' ,2gl5.7,//) 

Return to force subroutine. 

RETURN 

END 


ft ft ft ft ft ft' II ft ft * ft ft' ft ft ft ft ft ft ft ft ft ft. ft * * * ft * ft * f.-ft ft * ft ft ft ft mt ft- ft ft ft * fr.fl ftft ft 1* ftft .ft * ft ft ft ft ft ft ft ft ft ft ft 

SUBROUTI NE HSPACE ( ATRANS , ATEMP , ATEMP2 , ATEMP3 , BMATRX , 

1 C , CTEMP , FG , FPH I , FD , OUTPUT , C INVER , TEMPD , D ) 


Declare in COMMON all variables required by more than one 
subroutine. These variables are divided into four catagories: 

1. SIZE This common block contains variables which define 

array dimensions. 

2. PROB This common block contains variables which define 

problem parameters. 

3. MAP This common block contains variables which define 

. the finite element mesh. 

4. dbug This common block contains variables which define 

which debug switches are desired. 

Arrays are not passed through COMMON but passed as arguments to the 
subroutines, 


1 

1 


1 

1 

1 

1 


COMMON /SIZE/ NIJMANG, NUMRAD, NUMELM,INTNOD,NUMNOD,CODIAG, 

NUMEQN , LENGTH , WI DTH , SPACE , BNDELM , BCOLS , ALENG , AWIDE, ATMPL 
INTEGER NUMANG , NUMRAD , NUMELM , INTN0D , NUMNOD , COD I AG , NUMEQN , 

LENGTH , WIDTH , SPACE , BNDELM , BCOLS , ALENG , AWIDE , ATMPL 
COMMON /PROB/ W2 , GA , GB , RO , EXTG , ALPHA , THETA0 
REAL W2 , GA , GB , RO , EXTG , ALPHA , THETA0 

COMMON /MAP/ SWEEP, RANGE, DRAD, DANG, BDANG1 , BDANG2 , RAT IO, GRADNT 
INTEGER RATIO 

REAL SWEEP , RANGE , DRAD , DANG , BDANG1 , BDANG2 

LOGICAL GRADNT 

common /dbug/ dmain,dmesh,dassm,dktest,dforce, 

dhspac , drect , dtr i , dgamma , df ilam , damat , dbmat , 
dcmat , dnorm , ddiag , ddiag2 , df unc , dderiv , detest , datest 
logical dmain , dmesh , das sm , dktest , df orce , dh spac , 

drect , dtri , dgamma , df ilam , damat , dbmat , dcmat , dnorm , 
ddiag ,ddiag2 , df unc , dderiv , detest , datest 


Declare and dimension arrays used in this subroutine. 


OUTPUT ( I NTNOD) , BMATRX ( BNDELM , BCOLS ) 


REAL 
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COMPLEX ATEMP(ATMPL,*WIDE) ,CINVER(BNDELM, BNDELM) , 

1 TEKPD ( INTNOD, INTNOD) , D ( INTNOD, INTNOD) , 

1 ATRANS ( INTNOD, BNDELM) ,ATEMP2 (ALENG, AWIDE) , 

1 ATEMP3 (ALENG, AWIDE) ,FD( BNDELM) , 

1 C (BNDELM , 8NDELM) ,CTEMP(NUMANO) ,FG(INTNOD) ,PPHI (BNDELM) 

Declare variables used only in this subroutine. 

INTEGER R0W,C0L,I,PUSH,IER,C0L1,C0L2 
Determine B matrix. 

CALL BMAT(BMATRX, OUTPUT) 

Determine A matrix. 

CALL AMAT ( ATRANS , ATEMP , ATEMP2 , ATEMP3 ) 

If a test of the amatrix is desired run ATEST. 

IF (DATEST) CALL ATEST (ATRANS , FG , BMATRX) 

Complete A matrix by addition of term similar to B matrix. 


Iterate trough rows and columns of BMATRX, modify constant and add 
to appropriate term in ATRANS. 

DO 110 R0W=1, BNDELM 

DO 115 C0L=1,BC0LS 

PUSH=RATI0* (ROW-1) +C0L 
ATRANS (PUSH , ROW) = ATRANS (PUSH , ROW) - 
1 ALPHA* BMATRX ( ROW , COL ) 

115 CONTINUE 

110 CONTINUE 

++++++++++++++++++++++++++++++++++++++++++++++ 

If debug desired output ATRANS. 

if (.not.dhspac)goto 113 

write(26,5000) 

do ill row=l,intnod 

write (26,1000) (atrans(row,coi) ,col=l,bndeim) 

111 continue 

++++++++++++++++++++++++++++++++++++++++++++++ 

Determine C matrix. 

113 CALL CMAT (C , CTEMP ) 

If test of c matrix is desired run CTEST. 

I F ( DCTEST ) CALL CTEST (FPHI ,C) 

++++++++++++++++++++++++++++++++++++i+++++++++ 

If debug desired output C matrix. 

if (. not. dhspac) goto 36 
write (26,2000) 
do 35 row=l,bndelm 

write(26,l000) (c(row,col) ,col=l,bndelm) 

35 continue 
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C ++++++44+++++++++++++++++++++++++++++++++44-M-+ 

C 

c Perform condensation of exterior matrices to one matrix 
C representing the exterior at the boundary elements only, 
c 
c 

C Store C matrix in temporary storage so it can be inverted without 
C destruction. 

C 

36 DO 5 ROW"!, BNDELM 

DO 10 COL-1 , BNDELM 

ATEMP2(R0W,C0L)kC(R0W,CQL) 

10 CONTINUE 

5 CONTINUE 

invert c matrix according to instuctions in IMSL manual for LEQT1C, 

CALL LEQT1C (ATEMP2 , BNDELM , ALENG , CTEMP , 1 , BNDELM , 1 , OUTPUT , IER) 
CTEMP(l) = CMPLX(l. 0,0.0) 

DO 15 ROW=2, BNDELM 

CTEMP (ROW) - CMPLX(0. 0,0.0) 

15 CONTINUE 

DO 20 COL=l, BNDELM 

CALL LEQT1C ( ATEMP2 , BNDELM , ALENG , CTEMP , 1 , BNDELM , 2 , OUTPUT , I ER) 
DO 25 R0W=1, BNDELM 

CINVER(ROW,COL) - CTEMP (ROW) 

CTEMP (ROW) -• CMPLX(0. 0,0.0) 

IF ( ROW. EQ.COL+1) CTEMP (ROW) * CMPLX(l, 0,0.0) 

25 CONTINUE 

20 CONTINUE 

+++++++++++++++++++++++++++++++++++++++++++++++++++ 

If debug required output inverse of C matrix. 

if ( . not. dhspac) goto 31 
write (26,2500) 
do 30 row=l,bndelm 

write(26,1000) (cinver (row, col) ,col=l,bndelm) 

30 continue 

C +++++++++++++4"f++++++++++++++++++++++++++++++++++'f+ 

C 

C Initialize product of C inverse transpose times B to zero. 

C 

31 DO 40 R0W=1, BNDELM 

DO 45 COL=l , INTNOD 

TEMPD(ROW,COL) = CMPLX (0 . 0 , 0 . 0) 

45 CONTINUE 

40 CONTINUE 
C 

C Multiply C inverse transpose times B. Actual method is to multiply B 
C transpose times C inverse and then transposing the result. This is 
C done because of the way B is Stored. 

C 

DO 50 COLl=l, BNDELM 

DO 55 ROW=.l, BNDELM 

PUSH = RATIO* (ROW-1) 

DO 60 COL2=l,BCOLS 

TEMPD (COL1 , PUSH+C0L2 ) = 

1 TEMPD (COL1 , PUSH+COL2 ) + 

1 BMATRX ( ROW , COL2 ) * C I N VER ( ROW , COL1 ) 

60 


CONTINUE 
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CONTINUE 


55 

50 CONTINUE 

++++4 - +++++’t+++++++'t+'f++'t't r+’f++++4"f++++'f++4’‘f +++++■♦■ + 

If debug desired output C inverse transpose times B. 

if (. not. dhspac) goto 66 
write (26,3000) 
do 65 row*l, bndelm 

write (26,1000) (tempd(row,col) ,c,/,vi*l, intnod) 

65 continue 

+++ , f++++++4-++'f++++++++ - f+-f++++++++++++'f++'f4-4 , +++++++ 

Multiply A transpose times c inverse transpose times B. 

66 DO 70 CQL1 S 1, INTNOD 
DO 75 R0W«1, INTNOD 

ATEMP2 (ROW, COL1) * CMPLX(0. 0,0.0) 

DO 00 COL2*l, BNDELM 

ATEMP2 (ROW ,COLl) =ATEMP2 (ROW,COL1) + 

1 ATRANS ( ROW , COL2 ) * TEMPD (COL2 , COL1 ) 

80 CONTINUE 

75 CONTINUE 

70 CONTINUE 

+++++++++++++++++++++-+++++++++++++++■+++++++++++++++ 

If debug desired output A transpose times c inverse times B. 

if (.not.hspac) goto 86 
write (26, 3500) 
do 85 row=l, intnod 

write(26,1000) (ATEMP2 (row, col) ,col=l, intnod) 

85 continue 

+++ - f-+++++++'f+++++++++++')‘++++++++++++‘f ++++++++++++++ 

Complete half space matrix D by taking negative of the sum of the 
matrix just computed and its transpose. 

86 DO 90 R0W=1, INTNOD 
DO 95 COL*!, INTNOD 

D(R0W,COL)=-ATEMP2 (ROW, COL) -ATEMP2 (COL, ROW) 

95 CONTINUE 

90 CONTINUE 

+++++'f+-f+++++++++++++++++++++++++++++++-(-++++++ 

If debug desired output D matrix, 

if (.not. dhspac) goto 101 
write (26,4000) 
do 100 row-1, intnod 

write(26,1000) (d(row,col) ,col=l, intnod) 

100 continue 

C ++++++++ t+++++++++++++++++++++++++++++++++++++ 

c 

C Condense external force vectors to a single force vector at the 
c boundary, 
c 
c 

C Determine FD. 

C 

C Multiply B transpose times C inverse times FPHI. 

C 

101 DO 165 ROW=l, INTNOD 
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FD (ROW) - CMPLX(0. 0,0.0) 

DO 170 C0L*1,BNDELM 

FL'(ROW) ■ FD ( ROW ) +TEMPD ( COL , ROW ) » FPH I ( COL) 

170 CONTINUE 

165 CONTINUE 

c$ write(25,*) (fd(row) ,row*l,intnod) 
c 

C Multiply a transpose times C inverse transpose, 

c 

DO 175 C0L1*1,BNDELM 

DO 180 JROW*1,INTNOD 

TEMPD(ROW,COLl) ■ CMPLXfO. 0,0.0) 

DO 185 COL2«l,BNDELM 

TEMPD (ROW , COL1 ) « TEMPD(RQW,COLl) + 

1 ATRANS (ROW , COL2 ) «CI N VER ( COL1 , COL2 ) 

185 CONTINUE 

180 CONTINUE 

175 CONTINUE 

if ( .not .dhspac) goto 177 
write (26, 5500) 
do 176 row«i,intnod 

write(26,l000) (tempd(row,col) ,col“l,bndelm) 

176 continue 

c 

C Multiply A transpose times C inverse transpose times FPHI. 

C 

177 DO 190 ROW=l, INTNOD 

DO 195 COL-1, BNDELM 

FD(ROW) b FD (ROW) -TEMPO (ROW , COL) «FPHI (COL) 

195 CONTINUE 

190 CONTINUE 

c$ write(26,*) (fd(row) ^rowel.intnod) 
c 

C Debug format statements. 


c 

5500 FORMAT (/2X, **»***#»»* A TRANSPOSE TIMES C INVERSE TRANSPOSE 
1 

5000 FORMAT (/2X , '•»****»** ATRANS *»»******') 

4000 FORMAT (/2X, '***•»»*** D MATRIX **»*»»**«>) 

3500 FORMAT (/2X, * ««»****** ATRANS TIMES CINVER TIMES B *****•*•*') 
3000 FORMAT(/2X, ' ****»«**« C INVERS TRANSPOSE TIMES B «»**»»*«*•) 
2500 FORMAT (/2X, ’ *«*•*«*** C INVERSE »•«*»****') 

2000 FORMAT (/2X , 1 * * * *tm * * « C MATRIX *********' ) 

1000 FORMAT ( IX , 100F7 . 4 ) 

return 
END 

C »*»***»* ********************* *'******•* *:«*»**•* ********** 

C 

SUBROUTINE BMAT(BMATRX, OUTPUT) 

C 


C Declare in COMMON all variables required by more than one 
C subroutine. These variables are divided into four catagories: 

C 1. SIZE This common block contains variables which define 

C array dimensions. 

C 2. PROB This common block contains variables which define 

C problem parameters. 

c 3. MAP This common block contains variables which define 

C the finite element mesh. 

C 4. dbug This common block contains variables which define 
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which debug switches are desired. 

Arrays are not passed through COMMON but passed as arguments to the 
subroutines. 

COMMON /SIZE/ NUMANG,NUMRAD,NUMELM,INTNOD,NUMNOD,CODIAG, 

1 NUMEQN , LENGTH# WIDTH , SPACE , BNDELM , BCOLS , ALENG * AW1 DE * ATMPL 

INTEGER • NUMANG , NUMRAD , NUMELM , INTNOD , NUMHOD , COD I AG, NUMEQN , 

1 LENGTH .WIDTH# SPACE# BNDELM * BCOLS , ALENG , AWI DE , ATMPL 

COMMON /PROB/ W2 , GA , GB , RO , EXTG , ALPHA , THETAO 
REAL W2 , Qhi GB , RO , EXTG , ALPHA # THETAO 

COMMON /MV / SWEEP , RANGE * DRAD , DANG * BDANG.l , BDANG2 * RATIO * GRADNT 
INTEGER RATIO 

REAL SWEEP * RANGE * DRAD * DANG * BDANG1 # BDANG2 

LOGICAL GRADNT 

common /dbug/ dmain#dmesh,dassm#d)itest,dforce, 

1 dhspac , drect , d t r i , dgamma , df ilam , damat , dbmat , 

1 dcmat * dnorm # ddiag * ddiag2 ,df unc * dderiv .detest #datest 

logical dmain*dmesh,dassm#d)ctest,dforce, dhspac# 

1 drect , dt ri , dgamma , d f ilam , damat , dbmat , dcmat , dnorm * 

1 ddiag , ddiag2 * df unc # dderiv # detest , datest 

Declare and dimension arrays used in this subroutine. 

REAL BMATRX ( BNDELM , BCOLS ) , OUTPUT ( I NTNOD) 

Declare variables used in this subroutine only. 

INTEGER I, J, ELM# ROW, COL 
REAL ANS 

Calculate integral from derived expression. The terms of the integral 
for both shape functions are the same. 

ANS-0.25*DANG 
C 

c Position these values in the proper location of the B matrix. The 
C positioning is derived from the use of a circulant matrix which is 
C condensed to the one for a half-circle, only the non-zero terms are 
C stored. 

C 

DO 60 ROW=l, BNDELM 

BMATRX ( ROW , 1 ) =2 . 0 * ANS 
BMATRX ( ROW , BCOLS ) =2 . 0 * ANS 
IF (BCOLS. EQ. 2) GOTO 50 
DO 70 COL=2, BCOLS- 1 

BMATRX (ROW, COL) =2. 0*ANS 

70 CONTINUE 

60 CONTINUE 

C ++!+•+++++++++++++++++++++■♦-++++++++++++++++++ 

C If debug is desired output both the stored version of the B matrix and 
C the full representation for the half-circle. Otherwise return to 
c hspace subroutine, 
c 

if (.not. dbmat) RETURN 
write (26, 300) 
do 40 i=l,bndelm 

do 50 j=l#bcols 

output (j)*BMATRX(i, j) 

50 continue 

write(26,200) (output (j ) , j=l,bc<;is) 
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40 continue 

write (26,400) (i, i=l,intnod) 
do 10 row=l,bndelm 

do 20 col=l,intnod 

output(col)=0.0 

20 continue 

push * ratio* (row-1) 
do 30 00 . 1 = 1 , ratio+1 

output (push+col) =bmatrx ( row , col ) 

30 continue 

write (26, 100) row, (output (i) ,i=l,intnod) 

10 continue 

C +++++4 , +++++ , !'+++++++++++++++++++'H'++++++++++ 

C 

C Debug format statements 
C 

100 f ormat (2x, i2 , lx, lOOf 11 . 6) 

200 format(2x,10f6.3) 

300 fon;iat(/2y t , * ********* compact b matrix ****»*••••,/) 
400 format(//2x, ' ********* full b matrix **•»*••«**',// 
1 , 7:t, 100 (i2 , 4x) ) 

Return to HSPACE subroutine. 


return 

END 


a**************************************************************************** 


************************************ 

SUBROUTINE AMAT ( ATRANS , ATEMP , ATEMP2 , ATEMP3 ) 
C 


by more than one 
into four catagories: 
variables which define 


C Declare in COMMON all variables required 
C subroutine. These variables are divided 
C 1. SIZE This common block contains 

c array dimensions. 

C 2. PR0B This common block contains variables which define 

C problem parameters. 

C 3. MAP This common block contains variables which define 

c the finite element mesh. 

C 4. dbug This common block contains variables which define 

C which debug switches are desired. 

C Arrays are not passed through COMMON but passed as arguments to the 
C subroutines, 
c 


COMMON /SIZE/ NUMANG,NUMRAD, NUMELM , INTN0D, NUMNOD,CODIAG, 

1 NUMEQN , LENGTH , WIDTH , SPACE , BNDELM , BCOLS , ALENG , AWIDE , ATMPL 

INTEGER NUMANG , NUMRAD , NUMELM , INTNOD , NUMNQD , CODIAG, NUMEQN , 

1 LENGTH ,WIDTH, SPACE , BNDELM , BCOLS , ALENG , AWIDE , ATMPL 

COMMON /PROB/ W2 , GA , GB , RO , EXTG , ALPHA , THETA0 
REAL W2,GA,GB,RO, EXTG, ALPHA, THETA0 

COMMON /MAP/ SWEEP, RANGE ,DRAD, DANG, BDANG1 , BDANG2 , RAT IO, GRADNT 
INTEGER RATIO 

REAL SWEEP , RANGE , DRAD , DANG , BDANG1 , BDANG2 
LOGICAL GRADNT 

common /dbug/ dmain,dmeah,dassm,dktest,dforce, 

1 dhspac , drect , dtri , dgamma , df ilam , damat , dbmat , 

1 dcmat,dnorm,ddiag,ddiag2,dfunc,dderiv, detest, datest 

logical dmain , dmesh , dassm , dktest , df orce , dhspac , 

1 drect , dt r i , dgamma , df ilam , damat , dbmat , demat , dnorm , 

1 ddiag,ddiag2,dfunc,dderiv, detest, datest 


0 0-0 000 ooooooo ooo ooo ooo ooo 


Declare in COMMON the Gauss integration points. 

COMMON /INGRAL/ NUMPTY, YVALUE,WGTY ,NUMPTZ,ZVALUE,WGTZ 
INTEGER NUMPTY , NMPTYD , NUMPTZ 

Declare and dimension arrays required in this subroutine. 

COMPLEX ATEMP (2 , ALENG) ,TEMP, 

1 ATRANS ( I NTNOD , BNDELM } , ATEMP2 (ATMPL, AWIDE) ,ATEMP3 (ALENG, AWIDE) 
REAL YVALUE (100) ,WGTY (100) ,ZVALUE(100) ,WGTZ(1C0) 

Declare variables used only in this subroutine. 

INTEGER I , ROW , COL , SOURCE , PULL1 , PULL2 , OBSERV , PUSH , LIMIT , 

1 PULL , WOR , PUSH1 , PUSH2 

REAL QAOBS , QBOBS , QASRC , QBSRC , TEST 
LOGICAL ODD 

Read in Gauss integration points. 

READ (24,*) NUMPTY 
DO 5 I=l,NUMPTY/2 

READ (24, *) YVALUE (I) ,WGTY(I) 

5 CONTINUE 

READ (24, * ) NUMPTZ 
DO 10 I=l,FUMPTZ/2 

READ (24,*) ZVALUE ( I ) , WGTZ ( I ) 

10 CONTINUE 

+++++++++++++++++++++++++++++++++++++ 
if (damat ) write (26,7000) numptz , numpty 
++++++++++++++-f-4-+-f+++++++++++++++++++ 

Determine if number of interior boundary elements per exterior 
boundary element is odd or even. Set limits on the number of observer 
integrals necessary. 

ODD = .FALSE. 

LIMIT = RATIO/2 

TEST = RATIO - 2*LIMIT 

IF(TEST.NE.O) ODD = .TRUE. 

IF (LIMIT. EQ. 0) GOTO 16 

Perform double integral for A matrix. Outer integrals are the 
observers. 


Iterate over necessary observers. 

DO 30 0BSERV=1, LIMIT 
C +++++++++++++++++++++++++++++ 

if (damat)write(26,500) observ 
C +++++++++++++++++++++++++++++ 

c 

C Determine limits of integration for observer. 
C 

QAOBS = (OBSERV-1) *DANG 
QBOBS = OBSERV* DANG 

C +++++++++++++++++++-H-++++++++ 

if (damat ) write (26 , 600) qaObs ,qbobs 
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+++++++++++++++++++++++++++++ 

Iterate over necessary sources. 

DO 15 S0URCE=1 , BNDELM*2 

+++++++++++++++++++++++++++++ 
if (damat) write (26, 700) source 
+++++++++++++++++++++++++++++ 

Determine limits of integration for the source. 

QASRC = (SOURCE-1) *BDANG1 
QBSRC = S0URCE*BDANG1 
+++++++++++++++++++++++++++++ 
if (damat) write ( 26,800) qasrc,qbsrc 
+++++++++++++++++++++++++++++ 

Determine temporary storage slot for observer-source 
integral. 

PUSH = (OBSERV-l) *AWIDE+SOURCE 

Evaluate integral. 

IF (SOURCE. EQ.l) GOTO 11 

CALL N0RM2 (QAOBS , QBOBS , QASRC , QBSRC , 

1 ATEMP(1, PUSH) ,ATEMP(2, PUSH)) 

GOTO 12 

11 CALL NORM (QAOBS, QBOBS, QASRC, QBSRC, 

1 ATEMP (1 ,PUSH) ,ATEMP (2 ,PUSH) ) 

++++++++++++++++++++++++++++++++++++++++++++++ 

12 if (damat)write(26,900)push,atemp(l,push) ,push, 

1 atemp(2,push) 

++++++++++++++++++*+++++++++++++++++++++++++++ 

15 CONTINUE 
30 CONTINUE 

If number of internal boundary elements per external boundary element 
is odd perform integral for the odd observer. This is done here since 
the number of sources is less than for the other observers due to 
symmetry. 

IF( .NOT. ODD) GOTO 20 

16 OBSERV = LIMIT+1 

Determine limits of integration for observer. 

QAOBS = (OBSERV-l) *DANG 
QBOBS = OBSERV* DANG 
+++++++++++++++++++++++++++++++++ 
if (damat ) write (26 , 600) qaobs ,qbobs 
+++++++++++++++++++++++++++++++++ 

Iterate over necessary sources. 

DO 25 SOURCE=l , BNDELM+1 

+++++++++++++++-H-++++++++++++++++ 

if (damat) write (26, 700) source 

+++++++++++++++++++++++++++++-M-++ 
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Determine limits of integration for source. 

QASRC = (SOURCE-1) *BDANG1 
QBSRC = SOURCE* BDANG1 
+++++++++++ , f+4-++++++++++4' , f , f++++++ 

if (damat) write ( 26 , BOO) qasrc , qbsrc 

+++++++++++++++++++++++-f+++++++++ 

Determine temporary storage slot for observer-source 
integral. 

PUSH = (OBSERV-1) *AW1DE+S0URCE 
IF (SOURCE. EQ.l) GOTO 21 

Evaluate integral. 

CALL NORM2 (QAOBS , QBOBS , QASRC, QBSRC , 

ATEMP ( 1 , PUSH ) , ATEMP ( 2 , PUSH ) ) 

GOTO 22 

CALL NORM (QAOBS , QBOBS , QASRC , QBSRC, 

ATEMP (1, PUSH) , ATEMP (2, PUSH)) 
++++++++++++++++++++++++++++++++++++++++++++++ 
if (damat)write(26,900)push,atemp(l,push) ,push, 
atemp(2,push) 

+++++++++++++++++++++++++++++++++++•<•+++++++++4' 

CONTINUE 

4" 4- + 4- 4- +• 4" 4- + + 4- 4- 4" +■ 4" 4" 4-4 - 4- 4- 4- 4- 4-4- 4- 4- 4" 4- 4” 4- 4- 4- +■ 4- 4- 4- 4- 4- 4" 4- 4” 4" 4 , + 

If debug desired output computed integrals. 

20 if (.not. damat) goto 36 

write(26,300) 
limit=ratio»bndelm 

if (odd) limit=limit4-i 
do 35 col=l, limit 

write (26,400) (atemp(row,col) ,row=l,2) 

35 continue 

4-4‘4-4-4-4-4-+4-+4-4’4-4-4-4-+4-4-4-+4-4'4-4-4'4-4'+4-4-4'4-4-4-4-4-+4-4‘4 , 4-4-4- 

Determine if number of interior boundary elements per exterior 
boundary element is odd or even. Set limits on the length of A 
matrix and test value. 

36 LIMIT = RATIO/2 

IP (ODD) LIMIT = LIMIT+1 
TEST = RATIO *BNDELM 
IP (ODD) TEST=TEST4-1 


CREATE UNIT MATRIX . I 

DO 50 ROW*!, LIMIT j 

PUSH * ROW *2 

DO 40 COL=l,BNDELM*2 | 

PULL* (ROW-1) *BNDELM*2+COL - ' 

IF (PULL. GT. TEST) GOTO 37 ‘ 

ATEMP2 (PUSH-1 , COL) =ATEMP (1 , PULL) 

ATEMP2 (PUSH ,COL) =ATEMP (2 , PULL ) 

GOTO 40 . ' | 

37 PULL = 2*TEST-PULL j 

ATEMP2 ( PUSH-1 , COL) =ATEMP ( 2 , PULL) { 

ATEMP2 ( PUSH , COL ) =ATEMP ( 1 , PULL ) j 


j! 


l 

21 

1 

22 

1 

25 
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40 CONTINUE 

IF(ROW.EQ. LIMIT. AND. ODD)GOTO 50 

WOR = RATIO+l-ROW 

PUSH = W0R*2 

DO 45 COL=l , BNDELM* 2 

PULL=2+ROW* BNDELM* 2-COL 

I F ( PULL . GT . ROW * BNDELM * 2 ) PULL=PULL-BNDELM • 2 
ATEMP2 ( riJSH-1 , COL) =ATEMP ( 1 , PULL) 

ATEMP2 (PUSH , COL) =ATEMP (2 , PULL) 

45 CONTINUE 

50 CONTINUE 

if (.not.damat) goto 56 

write(26,2000) 

do 55 row=l,ratio*2 

write(26,7600) (atemp2 (row, col) ,col=l,bndelm*2) 

55 continue 

C SETUP BIG A MATRIX 

56 DO 70 1=1 , BNDELM*2 

DO 65 ROW=l , RATIO* 2 

PUSH1 = ROW+(I-l) *RATIO*2 
DO 60 COL=l,BNDELM*2 

PUSH2 = COL+I-1 

IF (PUSH2 .GT.BNDELM*2) PUSH2 = 

1 PUSH2-2* BNDELM 

ATEMP3 ( PUSH1 , PUSH2 ) =ATEMP2 (ROW , COL ) 

60 CONTINUE 

65 CONTINUE 

70 CONTINUE 

if ( .not.damat) goto 72 

write (26, 7100) 

do 71 row=l,bndelm*ratio«4 

write (26,7600) (atemp3 (row, col) ,col=l,bndelm*2) 

71 continue 

C ADD ROWS OF COMMON NODES 

72 DO 75 COL=l, BNDELM *2 

ATEMP3 (1 ,COL) = ATEMP3 (l,COL) +ATEMP3 (BNDELM*RATIO*4,COL) 

75 CONTINUE 

DO 85 ROW=2 , BNDELM*RATIO*4-2 
PULL = 2*ROW-2 
DO 80 COL=l, BNDELM*? 

ATEMP3 (ROW, COL) = ATEMP3 (PULL, COL) + 

1 ATEMP3 (PULL+l,COL) 

80 CONTINUE 

85 CONTINUE 

if (. not. ^T.'mat) goto 87 

write (26, 7200) 

do 86 row=l,bndelm*ratio*2 

write ( 26 , 7600 ) (atemp3 ( row , col ) , col=l , bndelm* 2 ) 

86 continue 

C FLIP RIGHT HAND GROUP OF COLUMNS 

87 DO 95 ROW=l,BNDELM*RATIO*2 

DO 90 1=1 , BNDELM/2 

PUSH1 = BNDELM+I 
PUSH2 = 2* BNDELM- ( I-l) 

TEMP = ATEMP3 (ROW , PUSH1 ) 

ATEMP3 (ROW , PUSH1) = ATEMP3 (ROW, PUSH2) 

ATEMP3 (ROW,PUSH2) = TEMP 

90 CONTINUE 

95 CONTINUE 

if (.not.damat) goto 97 


write (26, 7300) 

do 96 row=i, bndelm*ratio*2 

write(26,7600) (atemp3 (row, col) ,col=3,„bndelm*2) 

96 continue 

C FLIP LOWER GROUP OF ROWS 

97 DO 105 COL-1 , BNDELM * 2 

LIMIT = (BNDELM*RATI0)/2 
DO 100 1=1, LIMIT 

PUSH1 = BNDELM*RATIO+I 
PUSH2 = BNDELM*RATIO*2-(I-l) 

TEMP = ATEMP3 (PUSH1,C0L) 

ATEMP3 ( PUSH1 , COL) * ATEMP3(PUSH2,COL) 

ATEMP3 (PUSh'2,COL) = TEMP 

100 CONTINUE 

105 CONTINUE 

if ( . not. damat) goto 107 

write (26, 7400) 

do 106 row=l,bndelm»ratio*2 

write (26, 7600) (atemp3 (row, col) ,col=l,bndelm»2) 

106 continue 
C CONDENSE PHI 

107 DO 115 R0W=1 ,BNDELM*RATIO*2 

DO 110 COL=l, BNDELM 

PULL = BNDELM+COL 

ATEMP3 (ROW , COL) = ATEMP3 (ROW, COL) + 

1 ATEMP3 (ROW, PULL) 

110 CONTINUE 

115 CONTINUE 

if (.not, damat) goto 117 

write (26,7500) 

do 116 row=l,bndelm ,l| ratio*2 

write (26, 7600) (atemp3 (row, col) ,col=l,bndelm) 

116 continue 

C CONDENSE DISPLACEMENT 

117 DO 125 COL=l, BNDELM 

DO 120 ROW=2 , BNDELM* RATIO 

PULL = BNDELM*RATIO+ (ROW-1) 

ATEMP3 (ROW, COL) = ATEMP3 (ROW, COL) + 

1 ATEMP3 (PULL, COL) 

120 CONTINUE 

125 CONTINUE 

DO 130 COL=l, BNDELM 

ATEMP3 ( BNDELM * RAT 1 0 +1 ,COL) =ATEMP3 (BNDELM«RATIO*2 ,COL) 

130 CONTINUE 

DO 132 ROW=l, INTNOD 

DO 131 COL=l, BNDELM 

ATRANS ( ROW , COL ) =ATEMP3 ( ROW , COL ) 

131 CONTINUE 

132 CONTINUE 

if (.not. damat) RETURN 
write (26, 3000) 
do 13£j row=l,intnod 

write (26, 7600) (atrans (row, col) ,col=l,bndelm) 
if (row.eq.l.or .row. eq.bndelm* ratio) write (26, 6000) 

135 continue 

300 format (/,2x, '*****»» atemp transpose ******** * ,/) 

400 format (2x,100 (’(' ,2fll. 6, ') *)) 

500 format (2x, 'observ = • , ±5) 

600 format (2x,' qaobs=' ,fii.6, * qbobs-’ ,fll.6) 

700 format (2x,* source =',i5) 
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800 format(2x,' qasrc=' , fll.6, ’ qbsrc=* , fll.6) 

900 format(2x,' atemp(l, * ,13, ») =' # 2fll.6,/, 

1 2X, ' atemp(2, ' ,i3, ' ) =',2fll.6) 

1000 format(2X,8F6.3, ' |',8F6.3) 

2000 format(/,2x, ' *********** atemp2 *************) 

3000 format(/,2x, ' •*«•••«*•** atrans *•*****•«***•) 

5000 format(2x,8F6.3) 

6000 format (2x, ' ' ) 

7000 format (/,2x, 'number of points for dz integral =' ,i5,/, 

1 2x,'number of points for dy integral =',i5) 

7100 format(/,2x, ' ***«*»«* full a matrix **«****•') 

7200 format (/,2x, ' ******** adding common nodes ******** •) 

7300 format (/, 2x, ' ******** flipping righthand columns **»***«*') 

7400 format (/, 2x, ' ******** flipping lower rows ***•**«*») 

7500 format (/,2x, ' ******** condensing out phi *»**»***') 

7600 format (2x,20F6. 3) 
return 
END 

C********** **************** ^ ******************************************** 

c 

SUBROUTINE NORM (C , D , E, F , INTI , INT2) 

C 

C Declare in' COMMON all variables required by more than one 
C subroutine. These variables are divided into four catagoriess 
C 1. SIZE This common block contains variables which define 

c array dimensions. 

C 2. PROB This common block contains variables which define 

C problem parameters. 

C 3. MAP This common block contains variables which define 

c the finite element mesh. 

C 4. dbug This common block contains variables which define 

C which debug switches are desired. 

C Arrays are not passed through COMMON but passed as arguments to the 
c subroutines. 

I 

c 

COMMON /SIZE/ NUMANG , NUMRAD , NUMELM , I NTNOD , NUMNOD , COD I AG , 

1 NUMEQN , LENGTH /WIDTH , SPACE , BNDELM , BCOLS , ALENG , AWIDE , ATMPL 

INTEGER NUMANG , NUMRAD , NUMELM , INTNOD , NUMNOD , CODIAG , NUMEQN , 

1 LENGTH, WIDTH, SPACE, BNDELM, BCOLS, ALENG f AWIDE, ATMPL 

COMMON /PROB/ W2 , GA , GB , RO , EXTG , ALPHA , THETA0 
REAL W2 , GA , GB , RO , EXTG , ALPHA , THETA0 

COMMON /MAP/ SWEEP, RANGE, DRAD, DANG, BDANG1 , BDANG2 , RATIO, GRADNT 
INTEGER RATIO 

REAL SWEEP , RANGE , DRAD , DANG, BDANG1 , BDANG2 

LOGICAL GRADNT 

common /dbug/ dmain,dmesh,dassm,dktest,dforce, 

1 dhspac , drect , dtri , dgamma , df ilam , damat , dbmat , 

1 dcmat , dnorm , ddiag , ddiag2 , df unc » dderiv , detest , datest 

logical dmain ,dmesh ,dassm ,dktest ,dforce , dhspac , 

1 drect , dt ri , dgamma , df ilam , damat , dbmat , dcmat , dnorm , 

1 ddiag ,ddiag2 ,df unc , dderiv , detest , datest 

Declare in COMMON the Gauss integration points. 

COMMON /INGRAL/ NUMPT Y , Y VALUE , WGTY , NUMPTZ , ZVALUE , WGTZ 
INTEGER NUMPT Y,NMPT YD, NUMPTZ 

REAL YVALUE(IOO) ,WGTY (100) , ZVALUE (100) , WGTZ (100) 

Declare variables used only in this subroutine. 


OOO O O O O O O OOO O O O O OOO 


C 

INTEGER I,ZI,YI 

REAL A f B , C , D , Z , Y , EPSLON , DI ST , CONST 

COMPLEX INTGYZ, INTGY , ANSY , INTI , INT2 , LN ,X2 , DERIV 
LOGICAL DONE 

+++++++++++++++++++++++++++++++++++ 
if debug required output number of integration points. 

if (dnorm) write ( 26 , 100) numpty , numptz 
++++++++++++++++++++++++++++++++++ , t 

Define required constants. 

CONST = 1.0/ (2.0* (D-C)) 

Initialize integrals. 

INTI = CMPLX(0. 0,0.0) 

INT2 = CMPLX(0.0,0.0) 

Integrate with respect to dz over each element. 


Iterate through each Gauss integration point. 

DO 50 ZI=1 , NUMPTZ/2 

Transform -1 to 1 positive integration points to c to d points 

Z = (ZVALUE(ZI) * (D-C)+D+C)/2.0 
DONE = .FALSE. 

GOTO 30 

C Transform -1 to 1 integration points to c to d points. 

C This is done after the positive point has been integrated 

C with respect to dy. 

c 

25 Z= (-ZVALUE (ZI ) * (D-C) +D+C) /2 . 0 

DONE = .TRUE. 

C 

C For each value from the dz, integral integrate with respect 

C to dy. 

c 

C Initialize dy integral to zero. 

C 

30 INTGY = CMPLX(0. 0,0.0) 

C +++++++++++++++»! +.<•+++++++++ 

if (dnorm) write ( 26 , 500) zi,z 
C ++++++++++++++++{-++++++++++ 

C 

C Iterate through each Gauss integration point, 

c 

DO 40 Y 1=1 , N •gpTY/2 
C 

C The dy inte.tal is divided into two integrals to 

C handle the Ivgrithmic singularity. The first is from 

C the lower Limit of the source to the z value. 

C 

A - E 

a=z 

C Transient! -1 to 1 positive integration points to e to 
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z points. 

Y m (yVALUE(YI) » (B-A)+B+A)/2.0 
++++++++++++++++++++++++++++ 
if (dnorm)write(26,550) yi,y 
+++ 4 ++++++++++++++++++++++++ 

.Evaluate deriVitive of Green's Function at integration 
•points.. Multiply by scaling factor and sum for total 
dy integral t 

DIST = ABS (SIN ( (Z-Y)/2 .0) ) 

MSY ~ DERIV (DIST) * (B-A)/2 . 0 
•++++++++++++++++++++++++++++ 
if (dnorm) write ( 26 , 600) ansy 

4 +++++++++++++++++++++++++++ 

Transform' -1 to 1 negative integration points to e to 
z points.* 

t - (“YyALUE( Yl) * (B-A) +B+A)/2 . 0 
+• -f 4" + 4* 4" 4* + *f 4* + 4* +" + 4* 4* 4* 
if (dnorm) write (26,550) yi,y 
4 + 444 + 44 ++++++++++++++++++++ 

Evaluate derivative of Green's Function at integration 
points. Multiply by scaling factor and sum for total 
dv integral. 

DIST = ABS (SIN ( (Z-Y) /2. . 0) ) 

ANSY = ANSY + DERIV(DIST) * (B-A)/2.0 
++++++++++++++++++++++++++++ 
if (dnorm) write (26, 600) ansy 
++++++++++++++++++++++++++++ 

The second portion of the dy integral is from z to f. 

fc=Z 
B = F 

++++++++++++++++++++++++++++ 
if (dnorm) write (26, 550) yi,y 
++++++++++++++++++++++++++++ 

Transform -1 to 1 positive integration points to z to 
f points. 

Y a fXVALUE(Yl) *(B-A)+B+A)/2.0 

Evaluate derivitive of Green's Function at integration 
points. Multiply by scaling factor and sum for total 
dy integral. 

DIST = ABS (SIN ( (Z-Y)/2.0) ) 

ANSY = ANSY + DERIV (DIST) * (B-A) /2 . 0 
+++++++++++++++++++++++++++ 
if (dnorm)write (26,600)ansy 
+++++++++++++++++++++++++++ 

Transform -1 to 1 negative integration points to z to 
f points. 

Y = (-YVALUE(YI)*(B-A)+B+A)/2.0 
+++++++++++++++++++++++++++ 

if (dnorm)wbite(26,550) yi,y 
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++++++++++4"f++4"f+++'f++++++4 - 

Evaluate derivitive of Green's Function at integration 
points. Multiply by scaling factor and sum for total 
dy integral. 

DIET = ABS(SIN((Z-Y)/2.0)) 

ANSY *= ANSY + DERIV(DIST) * (B-A)/2.0 
+++++++++++++++++-H-++++++++ 

if (dnorm) write (26 , 600 ) ansy 
+++++++4-++++++++++ , f ++++++++ 

Multiply dy integral by the integration weights. 

INTGY *= INTGY + WGTY(YI) *ANSY 
++++++-f++++++++'t-+++++++‘f , t++ 
if (dnorm) write (26, 650) intgy 
+++++++++4-++++++^ : - s *•++++++ 

40 CONTINUE 

Multiply dy evaluation by dz weights and scaling factors 
and the required constants and sum to total dz integrals. 

INTGYZ = CQNST*WGTZ(Zl) * (D-C)/2.0*INTGY 
INTI = INTI + INTGYZ* (D-Z) 

INT2 = INT2 + INTGYZ* (Z-C) 

IF (DONE) GOTO 50 
GOTO 25 

50 CONTINUE 

+4-+++++++++++++++++++++++++++ 

If debug required output value of integral, 
if (dnorm) write (26,200) intgyz 
++++++++++‘t+++++++++ , f ++++++++ 

Debug format statements. 

100 format (/2x, 'number of points for DY integral =',i5,/, 

1 2x, 'number of points for DZ integral =' ,i5) 

200 format(2x, 'dz integral = ',2fll.6) 

500 format(2x, 'zi = ',i5,/,2x,'z = ’,fll.6) 

550 format (/,2x, 'yi = ',i5,/,2x,'y = »,fll.6) 

600 format (2x, 'ansy = ',2fll.6) 

650 format (2x, 'dy integral = ’,2fll.6) 

Return to AMAT subroutine. 

RETURN 

END 


SUBROUTINE NORM.? (C,D,E, F, INTI, INT2) 

C 

C Declare in COMMON all variables required by more than one 
C subroutine. These variables are divided into four catagories: 

C 1, SIZE This common block contains variables which define 

C array dimensions. 

C 2, PROB This common block contains variables which define 

C problem parameters. 

C 3. MAP This common block contains variables which define 

C the finite element mesh. 
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4. dbug This common block contains variables which define 
which debug switches are desired. 

Arrays are not passed through COMMON but passed as arguments to the 
subroutines. 

COMMON /SIZE/ NUMANG , NUMRAD , NUMELM , INTNOD , NUMNOD , COD1 AO , 

1 NUMEfiN , LENGTH , WIDTH , SPACE , BNDELM , BCOLS 

INTEGER NUMANG , NUMRAD , NUMELM , INTNOD , NUMNOD , CODI AG , NUMEQN , 

1 LENGTH , WIDTH, SPACE , BNDELM , BCOLS 

COMMON /PROB/ W2,GA,GB,RO,EXTG, ALPHA, THETAO 
REAL W2 , GA , GB , RO , EXTG , ALPHA , THETAO 

COMMON /MAP/ SWEEP, RANGE, DRAD, DANG, BDANG1,BDANG2, RATIO, GRADNT 
INTEGER RATIO 

REAL SWEEP , RANGE , DRAD , DANG , BDANG1 , BDANG2 
LOGICAL GRADNT 

common /dbug/ dmain,dmesh,dassm,dktest,dforce, 

1 dhspac , drect , dtri , dgamma , df ilam , damat , dbmat , 

1 dcmat, dnorm, ddiag,ddiag2,dfunc,dderiv, detest, datest 

logical dmain , dmesh , dassm , dkte st , df orce , dhspac , 

1 drect , dt ri , dgamma , df ilam , damat , dbmat , demat , dnorm , 

1 ddiag,ddiag2,dfunc,dderiv, detest, datest 

Declare in COMMON the Gauss integration points. 

COMMON /INGRAL/ NUMPTY , Y VALUE, WGTY , NUMPTZ, ZVALUE,WGTZ 
INTEGER NUMPTY, NMPTYD, NUMPTZ 

REAL YVALUE (100) ,WGTY(lOO) ,ZVALUE(100) ,WGTZ(100) 

Declare variables used only in this subroutine. 

INTEGER I,ZI,YI 

REAL A,B»C,D,Z,Y,E,F ,DI ST, CONST 

COMPLEX INTGYZ , 1NTGY , ANSY , INTI , INT2 , LN ,X2 , DERI V 

LOGICAL DONE 

++++++++++++++++++++++++++++++++++++ 

If debug required output number of Gauss integration points. 

if (dnorm) write (26,loo)numpty,numptz 
•f *}• 4* *f 4* + + + + *f + *f + 4* *f 4* + 4* + 4* + *f *f 4* 4* *f + + ■f + + + 4* 

Define required constant. 

CONST - 1.0/ (2.0* (D-C) ) 

Compute double Green's function integral for A matrix. 


Initialize integrals. 

INTI = CMPLX(0. 0,0.0) 

INT2 = CMPLX(0. 0,0.0) 

Integrate with respect to dz over each element. 


Iterate through each Gauss integration point. 

DO 50 Z 1=1, NUMPTZ/ 2 

Transform -1 to l positive integration points to c to d points. 
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c 

Z * (ZVALUB(ZI) • (D-C)+D+C)/2 *0 
DONE « .FALSE. 

00X0 30 
C 

c Transform "1 to l negative integration points to c to d 

C points. This is done after the positive point has been 

C integrated with respect to dy. 

c 

25 Z=* (-ZVALUE(ZI) * (D-C)+D+C)/2.0 

DONE * .TRUE. 

C 

c For each value from the dz integral integrate with respect 

c to dy. 

c 
c 

c Initialize dy integral to zero, 

c 

30 INTGY = CMPLX(0.0;0.0) 

Iterate through each Gauss integration point. 

DO 40 YT=l,NUMPTY/2 

Define limits of integration. 

A « E 
B = F 

Transform -1 to 1 positive integration points to e to 
f points. 

Y = (YVALUE(YI) * (B-A)+B+A)/2.0 

Evaluate derivative of Green's function at integration 
points. Multiply by scaling factor and sum for total 
dy integral. 

DIST = ABS (SIN < (Z-Y)/2 .0) ) 

ANSY = DERIV (DIST) * (B-A)/2 .0 

Transform -1 to 1 negative integration points to e to 
f points. 

Y = (-YVALUE(YI) * (B-A) +B+AJ/2 . 0 

Evaluate derivative of Green's function at integration 
points. Multiply by scaling factor and sum for total 
dy integral. 

DIST = ABS (SIN ( (Z-Yl/2 .0) ) 

ANSY = ANSY + DERIV (DIST) * (B-A) /2 • 0 

Multiply dy integral by integration weights. 

INTGY = INTGY + WGTY(Yl) *ANSY 

40 CONTINUE 

Multiply dy evaluation by dz weights and scaling factors 
and the required constants and sum to tot&l dz integrals. 
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INTGYZ * CONST*WGTZ(ZI)*(D-C)/2.0*INTGY 
INTI * INTI + INTGYZ* (D-Z) 

INT2 * INT2 + INTGYZ* (Z-C) 

IF (DONE) GOTO 50 
GOTO 25 

50 CONTINUE 

++++++++++++++++++++++ 4 -++ 4 -++++++ 

If debug required output values of integrals. 

if (dnorm)vrite(26,200) intgyz 
if (dnorm) write ( 26 , 1100) inti , int2 
++++++++++++++++++++++++++++++++ 

Debug format statements. 

100 format(/2x, number of points for DY integral <=' ,i5,/, 

1 2x,' number of points for DZ integral =',15) 

200 format (2x, * integral = ' ,2f.ll.6) 

.100 format (2X, ' inti =' ,2fll,6,/, 

2x,' infc2 s',2fll.6) 

Return to AMAT subroutine. 


RETURN 

END 

******* •****$; 


* 

COMPLEX FUNCTION DERIV(DIST) 

C 


C Declare in COMMON all variables required by more than one 
C subroutine*. These variables are divided into four catagories: 
c 1. SIZE This common block contains variables which define 

C array dimensions. 

C 2. PROB This common block contains variables which define 

c problem parameters . 

C 3. MAP This common block contains variables which define 

C the finite element mesh. 

C 4. dbug This common block contains variables which define 

C which debug switches are desired. 

C Arrays are not passed through COMMON but passed as arguments to the 
C subroutines, 
c 


COMMON /SIZE/ NUMANG , NUMRAD , NUMELM , INTNOD , NUMNOD , CODIAG , 
1 NUMEQN , LENGTH .WIDTH , SPACE , BNDELM , BCOLS , ALENG , AW I DE , ATMPL 

I NTEGER NUMANG , NUMRAD , NUMELM , INTNOD , NUMNOD , CODIAG , NUMEQN , 

1 LENGTH , WIDTH , SPACE , BNDELM , BCOLS , ALENG , AWI DE , ATMPL 

COMMON /PROB/ W2,GA,GB,RO,EXTG, ALPHA, THETA0 
REAL W2 , GA , GB , RO , EXTG , ALPHA , THETA0 

COMMON /MAP/ SWEEP, RANGE, DRAD, DANG, BDANG1.BDANG2, RATIO, GRADNT 
INTEGER RATIO 

REAL SWEEP , RANGE , DRAD , DANG , BDANG1 , BDANG2 

LOGICAL GRADNT 


1 

1 

1 

1 


common /dbug/ dmain,dmesh,dassm,dktest,dforce, 

dhspac , drect , dtri , dgamma , df Ham , damat , dbmat , 
dcmat , dnorm, ddiag ,ddiag2 ,dfunc ,dderiv, detest ,datest 
logical dmain,dmesh,dassm,dktest,dforce, dhspac, 

drect , dt ri , dgamma , df i lam , damat , dbmat , dcmat , dnorm , 
ddiag, ddiag2,dfunc,dderiv, detest, datest 


non o ooo non nno o ooo ooo ooo ooo ooo ofio 


Declare variables used In this function only. 

REAL KAPPA , X , D1 ST , ARG , FI , THETA , Y 1 , J1 

Define necessary constant. 

KAPPA - SQRT(W2 GRANGE/ EXTG) 

Compute argument of Green's function derivative. 

X « 2 .0*KAPPA*RANGE*DIST 

Test argument to determine proper formula for computation. 

IP(X.LT.3.0) GOTO 10 
Compute derivative of Green's function. 

ARG * 3.0/X 

FI * ((((<-. 00020033»ARG+. 00113653) «ARG-. 00249511) *ARG 
1 +. 00017105) * ARG+. 01659667) *ARG+ . 00000156) *ARG 

1 +.79788456 

THETA « (((({-. 00029166* ARG+. 00079824) »ARG+. 00074348) *ARG 
1 - . 00637879) *ARG+. 00005650) • ARG+ . 12499612) • ARG 

1 -2.35619449+X 

J1 « FI "COS (THETA) /SgRT (X) * KAPPA* DIST/ (-4 . 0) 

Y1 “ F1*SIN (THETA) /SQRT(X) «KAPPA*DIST/ (4.0) 

+++ ++++++++++++++++++++++++++++++++++++++ 

If debug desired output computed values. 

i,f (dderiv) write (26, 500) kappa, dist,x, jl,yl 
+++++++++++++++++++++++++++++++++++++++++ 

DERIV = CMPLX( Yl, Jl) 

Return to NORM or N0RM2 subroutines. 

RETURN 

Compute derivative of Green's function. 

10 ARG = (X/3 .0) **2 

FI = ((((((. 00001109*ARG-. 00031761) *ARG+. 00443319) *ARG 
1 - . 03954289) *ARG+. 21093573) *ARG- . 56249985) «ARG+0. 5) *X 

THETA = ((((((. 0027873*ARG-. 0400976) *ARG+. 3123951) "ARG 
1 -1 . 3164827 ) » ARG+2 . 1682709) • ARG+ . 2212091 ) » ARG 

1 - . 6366198+2/3 . 1415927*X»ALOG (0 . 5»X) *F1) /X 

Jl = Fl*KAPPA«DIST/(-4.0) 

Y1 = THETA«KAPPA*DIST/(4.0) 
+++++++++++++++++++++++++++++++++++++++++ 

If debug desired output computed values. 

if (dderiv) write (26, 500) kappa, dist ,x, jl ,yl 
+++++++++++++++++++++++++++++++++++++++++ 

DERIV = CMPLX(Y1, Jl) 

Return to NORM or NORM 2 subroutines. 

RETURN 
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Debug format statements. 

500 format (2x, 'kappa = 

1 2x, 'distance = *,fll.6,/, 

1 2x, 'x » * ,fll.6,/, 

1 2X,'jl S ' , f 11.6,/ * 

1 2x, 'yl s »,fll.6,/, 

1 2x»'func = cmplx(yl, jl) ' ) 
END 


# 


SUBROUTINE ATEST (ATRANS , FG, BMATRX) 

C 

C Declare in COMMON all variables required by more than one 
c subroutine. These variables are divided into four catagories: 

C 1. SIZE This common block contains variables which define 

C array dimensions. 

C 2. PROB This common block contains variables which define 

C problem parameters . 

C 3. MAP This common block contains variables which define 

C the finite element mesh. 

C 4. dbug This common block contains variables which define 

C which debug switches are desired. 

C Arrays are not passed through COMMON but passed as arguments to the 
C subroutines. 


c 


1 

1 


1 

1 


COMMON /SIZE/ NUMANG , NUMRAD , NUMELM , INTNOD , NUMNOD , CODIAG , 
NUMEQN , LENGTH, WIDTH , SPACE , BNDELM , BCOLS , ALENG , AW IDE , ATMPL 
INTEGER NUMANG , NUMRAD , NUMELM , INTNOD , NUMNOD , CODIAG , NUMEQN , 

LENGTH , WIDTH , SPACE , BNDELM , BCOLS , ALENG , AWIDE , ATMPL 
COMMON /PROB/ W2 , GA , GB , RO , EXTG , ALPHA , THETAO 
REAL W2 , GA , GB , RO , E&TG , ALPHA , THETAO 

COMMON /MAP/ SWEEP, RANGE, DRAD, DANG, BDANG1 , BDANG2 , RATIO, GRADNT 
INTEGER RATIO 

REAL ' SWEEP, RANGE, DRAD, DANG, BDANG1 , BDANG2 
LOGICAL GRADNT 

common /dbug/ dmain,dmesh,dassm,dktest,dforce, 

dhspac , drect , dtri , dgamma , df ilam , damat , dbmat , 
dcmat, dnorm,ddiag,ddiag2,dfunc,dderiv, detest, datest 
logical dmain , dmesh , das sm , dktest , df orce , dhspac , 


1 drect , dtr i , dgamma , df ilam , damat , dbmat , dcmat , dnorm , 

1 ddiag,ddiag2,dfunc,dderiv, detest, datest 

INTEGER COL1 , COL2 , ROW , I ER 

REAL WA (15) , BMATRX (BNDELM, BCOLS) ,TMPB(15,15) ,TEMPB (15 , 16) 
COMPLEX ATRANS ( INTNOD, BNDELH) , 

1 FG( INTNOD) ,TMPA(15»15) ,TMPFG(15) ,ANS(15) ,TEMPA(15,16) 

1 ,TMPAT(16,15) ,TMPAAT(15,15) 

WRITE (26, 2000) 

I F ( I NTNOD . GT . 16 . OR . BNDELM . GT . 15 ) GOTO 100 
C Write BMATRX 

WRITE (26, 1100) 

DO 5 ROW=l, BNDELM 

WRITE(26,1000) (BMATRX(ROW,COLl) ,C0L1=1, BCOLS) 

5 CONTINUE 
C WRITE ATRANS 

WRITE (26, 1200) 

DO 10 ROW=l, INTNOD 

WRITE(26,1000) (ATRANS (ROW, COL1) , COL? =1, BNDELM) 

10 CONTINUE 


C COPY ATRANS 

DO 12 ROW-1, INTNOD 

DO 11 COLl=l , BNDELM 

TMPAT < ROW , COLl ) *ATRANS ( ROW , COLl ) 

11 CONTINUE 

12 CONTINUE 

C WRITE COPY OF ATRANS 
WRITE (26, 1900) 

DO 13 ROW=l, INTNOD 

WRITE (26, 1000) (TMPAT (ROW, COLl) ,COLlnl , BNDELM) 

13 CONTINUE 

C ADD BMATRX TERM TO ATRANS 
DO 20 ROW=l, BNDELM 

DO 15 COL1-1 * BCOLS 

COL2=RATIO* (ROW-1) +COL1 

TMPAT (COL2 , ROW ) =ATRANS (COL2 , ROW) + 

1 EXTG*ALPHA/G* BMATRX (ROW, COLl) 

15 CONTINUE 

20 CONTINUE 

C WRITE ATRANS WITH BMATRX TERM ADDED 
WRITE(26,.l300) 

DO 25 ROW=l, INTNOD 

WRITE (26, 1000) (TMPAT (ROW, COLl) , C0L1=1 , BNDELM ) 

25 CONTINUE 

C MULTIPLY A * A TRANSPOSE 
DO 40 C0L1=1, BNDELM 

DO 35 ROW=l, BNDELM 

TMPAAT (ROW , COLl) =CMPLX (0 .0,0.0) 

DO 30 COL2=l, INTNOD 

TMPAAT (ROW , COLl ) =TMPAAT (ROW , COLl ) + 

1 TMPAT (COL2 , ROW) * TMPAT (COL2 , COLl ) 

30 CONTINUE 

35 CONTINUE 

40 CONTINUE 
C WRITE A * A TRANSPOSE 
WRITE (26, 1400) 

DO 45 ROW=l, BNDELM 

WRITE (26, 1000) (TMPAAT (ROW, COLl) ,C0L1=1, BNDELM) 

45 CONTINUE 
C MULTIPLY A • FG 

DO 55 ROW=l, BNDELM 

TMPFG (ROW) =CMPLX ( 0 < 0,0.0) 

DO 50 COL2=l, INTNOD 

TMPFG (ROW) =TMPFG (ROW) + 

1 TMPAT (COL2, ROW) *FG(COL2) 

50 CONTINUE 

55 CONTINUE 

C WRITE A » FG 

WR T TE(26,1500) 

DO 60 ROW=l, BNDELM 
WRITE (26, 1000) TMPFG (ROW) 

60 CONTINUE 
C SOLVE FOR PHI 

CALL LEQT1C (TMPAAT , BNDELM , 15 , TMPFG , 1 , 15 , 0 , WA , IER) 

IF(IER.EQ.129) GOTO 110 
C WRITE PHI 

WRITE (26, 1600) 

DO 65 ROW=l, BNDELM 

WRITE (26, 1000) TMPFG (ROW) 

65 CONTINUE 
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VRITE(26,2000) 

RETURN 

100 WRITE(26,1700) 

STOP 

110 WRITE*26,180Q) 

STOP 

2000 FORMAT*/, 2X, ’ATEST SUBROUTINE ') 

1900 FORMAT*/, 2X, '*«**«*»** COPIED ATRANS *«***-«*.*') 

1800 FORMAT */,2X, 'A * A TRANSPOSE IS SINGULAR. ' ) 

1700 FORMAT */,2X, 'MATRICES ARE LARGER THAN THE TEMPORARY STORAGE 

1 IN ATEST. ' ,//,2X,' EXECUTION ABORTED') 

1600 FORMAT*/ ,2X, '********«* PHI «••**•***•) 

1500 FORMAT*/, 2X, '****•*»*»* A * FG *••*»**«*<.') 

1400 FORMAT */,2X, '««*»*i***m A*A TRANSPOSE ***«**«*»') 

1300 FORMAT*/, 2X , '***••*•«•* ATRANS + B TERM ****«**»*') 

1200 FORMAT */,2X, '»«****»*** ATRANS «***«***') 

1100 FORMAT */,2X, '***• ****** BMATRX ft*******') 

1000 FORMAT *2X, 100F11 . 6) 

END 


C 

C** *********************** K **************************************** * 


SUBROUTINE CMAT(C,CTEMF) 

C 

C Declare in COMMON all variables required by more than one 
c subroutine. These variables are divided into four catagories: 

C 1. SIZE This common block contains variables which define 

C array dimensions. 

C 2. PROB This common block contains variables which define 

C problem parameters. 

C 3. MAP This common block contains variables which define 

C the finite element mesh. 

c 4. dbug This common block contains variables which define 

C which debug switches are desired. 

C Arrays are not passed through COMMON but passed as arguments to the 
C subroutines. 


C 


1 


1 


1 


1 

1 


COMMON /SIZE/ NUMANG,NUMRAD,NUMELM , INTNOD,NUMNOD,CODIAG, 

NUMEQN , LENGTH , WIDTH , SPACE , BNDELM , BCOLS , ALENG , AW IDE , ATMPL 
I NTEGER NUMANG , NIJMRAD , NUMELM , I NTNOD , NUMNOD , COD I AG , NUMEQN , 

LENGTH , WIDTH , SPACE , BNDELM , BCOLS , ALENG , AWIDE , ATMPL 
COMMON /PROB/ W2 , GA , GB , RO , EXTG , ALPHA , THETAO 
REAL W2 , GA , GB , RO , EXTG , ALPHA , THETAO 

COMMON /MAP/ SWEEP, RANGE, DRAD, DANG, BDANG1 , BDANG2 , RAT 10, GRADNT 
INTEGER RATIO 

REAL SWEEP , RANGE , DRAD , DANG , BDANG1 , BDANG2 
LOGICAL GRADNT 

common /dbug/ dmain,dmesh,dassm,dktest,dforce, 

dhspac ,drect,dtri, dgamma , df ilarn , damat , dbmat , 
dcmat , dnorm , ddiag , ddiag2 ,df unc , dderiv , detest ,datest 
logical dmain,dmesh,dassm,dktest,dforce, dhspac, 

drect , dtri , dgamma , df ilarn , damat , dbmat , dema t , dnorm , 
ddiag , ddiag2 , df unc , dderiv , detest , datest 


Declare in COMMON the Gauss integration points. 

COMMON /INGRAL/ NUMPTY , Y VALUE , WGTY , NUMPTZ , ZVAL'JE , WGTZ 
INTEGER NUMPTY, NMPT YD, NUMPTZ 

REAL YVALUE (100) , WGTY (100) ,ZVALUE(100) ,WGTZ(100) 


ono ooooono ooo oo 


Declare and dimension arrays required in this subroutine. 

COMPLEX C ( BNDELM , BN DELM ) ,CTEMP (NUMANG) 

Declare variables used only in this subroutine. 

I NTEGER I , ROW , COL , SOURCE , PULL1 , PULL2 
REAL QAOBS , QBOBS , QASRC , QBSRC 

Perform double integral of Green's function. 


Determine limits of integration for outer integral. This is the 
integral over the observer. 

QAOBS =0.0 
QBOBS = BDANG1 

Perform integral over first source. 

SOURCE = 1 
C 

C Read in the Gauss integration points desired. 

C 

READ ( 24 , * ) NUMPTY 
DO 10 I=l,NUMPTY/2 

READ<24,») YVALUE(I) ,WGTY(I) 

10 CONTINUE 

READ (24 , * ) NUMPTZ 
DO 20 1=1 / NUMPTZ/2 

READ (24 , * ) Z VALUE ( I ) , WGTZ ( I ) 

20 CONTINUE 

C 

C Evaluate integral. 

C 

CALL D I AG ( QAOBS , QBOBS , CTEMP ( SOURCE ) ) 

C +++++++++++++++++++++++++++++++++-M-++++++++++++++++++++++++ 

C If debug desired compute exact integral for ln|z-y| and 

C output. 

C 

if (dcmat)exact=(qaobs-qbobs) * (qbobs-qaobs) + (qbobs-qaobs) **2 
1 * (alog(qbobs--qaobs)-0.5)+qbobs*(1.0-alog(qbobs-qaobs) ) * 

1 (qbobs-qaobs-1.0) 

if (dcmat) write (26,600) exact 
if (dcmat)write (26,500) source, ctemp (source) 

C ++++++++++++++++++++++++++++++++++++++++-«-+++++-t-++++++++++++ 

C 

C Perform integral over second source. 

C 

SOURCE = 2 
C 

C Read in the Gauss integration points desired. 

C 

READ (24,*) NUMPTY 
DO 40 I=l,NUMPTY/2 

READ(24, *) Y VALUE ( I ) , WGTY ( I ) 

40 CONTINUE 

READ (24,*) NUMPTZ 
DO 50 1=1, NUMPTZ/2 

READ (24,*) Z VALUE (I), WGTZ (I) 
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50 CONTINUE 

Determine limits of integration for second source. 

QASRC = ( SOURCE-1 )*BDANG1 
QBSRC = SOURCE* BDANG1 

Evaluate integral 

CALL DIAG2 (QAOBS , QBOBS , QASRC , QBSRC , CTEMP (SOURCE) ) 
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
If debug desired compute exact integral for ln|z-y| and 
output. 

if (dcmat) exact®-. 5* (qbsrc-qbobs) *«2*alog (qbsrc-qbobs)+. 75* 

1 (qbsrc-qbobs) **2+.5*qbsrc**2*alog (qbsrc)-. 75*qbsrc**2+ 

1 -1 . 0* ( . 75* 

1 (qasrc-qbobs) «*2+.5*qasrc* *2 *alog (qasrc) 73»qasrc**2) 
if (dcmat)write (26,600) exact 
if (dcmat ) write (26,500) source ,ctemp ( source) 
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 

Read in Gauss integration points for remaining sources. 

READ (24 , * ) NUMPTY 
DO 60 1=1 ,NUMPTY/2 

READ(24, *) YVALUE(I) ,WGTY(I) 

60 CONTINUE 

READ (24, * ) NUMPTZ 
DO 70 1=1 , NUMPTZ/2 

READ (24, *) ZVALUE ( I ) , WGTZ ( I ) 

70 CONTINUE 

Perform integral over remaining sources. 

DO 80 s'oURCE=3 , BNDELM + 1 

Determine limits of integration. 

QASRC = (SOURCE-1) *BDANG1 
QBSRC = SOURCE* BDANG1 

Evaluate integral. 

CALL DIAG2 (QAOBS, QBOBS, QASRC, QBSRC, CTEMP (SOURCE) ) 

C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 

C If debug desired compute exact integral for Injz-yj and 

C output. 

C 

if (dcmat)exact=-.5* (qbsrc-qbobs) * *2 *alog (qbsrc-qbobs )+. 75* 

1 (qbsrc-qbobs) **2+.5*qbsrc**2»alog (qbsrc)-. 75*qbsrc**2+ 

1 -1.0* (- .5* (qasrc-qbobs) **2*alog (qasrc-qbobs) +.75* 

1 (qasrc-qbobs) **2+.5*qasrc**2*alog(qasrc)-.?5*qasrc**2) 

if (dcmat)write (26,600) exact 
if (dcmat) write (26, 500) source, ctemp (source) 

C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 

80 CONTINUE 

C ++++-H+++++++++++++++++++++++++++++++ 

if (dcmat) write (26,400) (i,i=l,bndelm) 

c +++++++++++++++++++++++++++++++++++++ 
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C 

C Position observer one integrals into proper place in C matrix. This is 
c done on the basis of a circulant matrix which is then condensed to 
c reflect only a half-circle. 

C 

DO 100 R0W=1, BNDELM 

DO 90 C0L«=1, BNDELM 

IF(COL.LT.ROW) GOTO 05 
PULL1 = COL- (ROW-1) 

PULL2 = COL + ROW 

IF (PULL2 .GT . BNDELM+1) PULL2=2« (BNDELM+1) -PULL2 
C (ROW, COL) = 2 . 0» (CTEMP (PULL1) + CTEMP(PULL2) ) 
GOTO 90 

85 0 (ROW, COL) = C (COL, ROW) 

90 CONTINUE 

C ++++++++++++++++++++++++++++++++++++++++++++++++++++ 

if (dcmat)write (26, 300) row, (c(row,col) ,col=l,bndelm) 

C ' ++++++++++++++++++++++++++++++++++++++++++++++++++++ 

100 CONTINUE 

Debug format statements. 

300 fora-?,' , „X,i2,lX,100 (lX,2f 6. 3) ) 

400 format(//2x, '««***»*»* full c matrix «*••*«**•*' ,// 

1 , llx,100(i2,llx) ) 

500 format (2x, 'ctempC ,i3, ' ) = ' ,2fll.6,/) 

600 format (2x, 'exact ln|z-y| integral = ' ,fll.6) 

Return to HSPACE subroutine. 

RETURN 
END 


SUBROUTINE DlAG(C,D, INTGYZ) 

Declare in COMMON all variables required by more than one 
subroutine. These variables are divided into four catagories: 

1. SIZE This common block contains variables which define 

array dimensions. 

2. PROB This common block contains variables which define 

problem parameters. 

3. MAP This common block contains variables which define 

the finite element mesh. 

4. dbug This common block contains variables which define 

which debug switches are desired. 

Arrays are not passed through COMMON but passed e.s arguments to the 
subroutines. 

COMMON /SIZE/ NUMANG , NUMRAD , NUMELM , INTNOD , N'JMNOD , C0DIAG , 

1 NUMEQN , LENGTH , WIDTH , SPACE , BNDELM , BCOLS , ALENG, AWIDE, ATMPL 

INTEGER NUMANG , NUMRAD , NUMELM , INTNOD , NUMNOD , COD I AG, NUMEQN , 

1 LENGTH, WIDTH, SPACE, BNDELM, BCOLS, ALENG, AWIDE, ATMPL 

COMMON /PROB/ W2 , GA , GB , R0 , EXTG , ALPHA , THETA0 
REAL W2 , GA , GB , RO , EXTG , ALPHA , THETA0 

COMMON /MAP/ SWEEP, RANGE, DRAD, DANG, BDANG1 , BDANG2 , RATIO, GRADNT 
INTEGER RATIO 

REAL SWEEP , RANGE , DRAD , DANG , BDANG1 , BDANG2 

LOGICAL GRADNT 


ooo ooooan ooooao 0000 ooo ooo ooo 


1 

1 

1 

1 


common /dbug/ dmain,dmesh,dassm,d)ctest,dforce, 

dhspac , d r e c t , d t ri , dgamma , df il am , damat , dbmat , 
dcmat , dnorm , ddiag , ddiag2 , df unc , dderiv .detest, dates t 
logical dmain , dmesh , dassm, dktest , df orce , dhspac , 

drect.dtri, dgamma , df i lam , damat , dbmat , dcmat , dnorm , 
ddiag ,ddiag2 ,df unc , dderiv, detest ,datest 


Declare in COMMON the Gauss integration points. 

COMMON /INGRAL/ NUMPTY, YVALUE,WGl!Y,NUMPTZ,ZVALUE,WGTZ 
INTEGER NUMPTY »NMPTYD,NUMPTZ 

REAL YVALUE(IOO) ,WGTY(100) ,ZVALUE(100) ,WGTZ(100) 
Declare variable used only int this subroutine. 

INTEGER I,ZI,YI 

REAL A,B,C,D,Z,Y , EPSLON ,D 1ST, CONST 

COMPLEX INTGYZ, INTGY ,ANSY ,FUNC 
LOGICAL DONE 

++++++++^ (•+++++++++++++++++++++++++ 

If debug desired output number of Gauss integration points. 


if (ddiag) write (26,100) numpty , numptz 

++++++++++ H++++++++++++++++++++++++ 

Define required constant 
CONST = -0.5 

Compute double Green's function integral for C matrix. 


Initialize integral 

INTGYZ = CMPLX(0. 0,0.0) 


Integrate with respect to dz over each element. 
Iterate through each Gauss integration point. 


DO 50 ZI=1 , NUMPTZ/2 

Transform -1 to 1 positive integration points to c to d points 


c 

c 

c 

c 

c 


25 


C 

C 

C 

C 

C 


Z = ( Z VALUE (ZI ) * (D-C) +D+C) /2 . 0 
DONE “ .FALSE. 

GOTO 30 

Transform -1 to 1 negative integration points to c to d 
points. This is done after the positive point has been 
integrated with respect to dy. 

Z= (-ZVALUE (ZI ) * (D-C) +D+C) /2 . 0 
DONE = .TRUE. 

For each value from the dz integral integrate with respect to dy. 


Initialize dy integral to zero. 


ooo ooooo oooo ooooo o o o o ooooo o o o 


C 

30 INTGY = CMPLX(0. 0,0.0) 

Iterate through each Gauss integration point. 

DO 40 YI=l,NUMPTY/2 

The dy integral is divide into two integrals to handle 
the logrithmic singularity. The first is from the lower 
limit to the z value. 


A = C 
B=Z 

Transform -1 to 1 positive integration points to c to 
z points. 

Y = (YVALUE(YI) * (B~A)+B+A)/2.0 

Evaluate Green's function at integration points. 
Multiply by scaling factor and sum for total dy 
integral. 

DIST = 2.0«RANGE*ABS (SIN ( (Z-Y)/2 .0) ) 

ANSY = FUNC(DIST)* (B-A)/2.0 

Transform -1 to 1 negative integration points to c to z 
points. 

Y = (-YVAL'JE(YI) » (B-A)+B+A)/2.0 
DIST = 2 . 0* RANGE* ABS (SIN ( (Z-Y) /2 . 0) ) 

Evaluate Green's function at integration points. 
Multiply by scaling factor and sum for total dy 
integral. 

ANSY = ANSY + FUNC (DIST) * (B-A) /2 . 0 

The second portion of the dy integral is from z to d. 

A=Z 
B = D 
C 

C Transform -1 to 1 positive integration points to c to 

C z points. 

C 

Y = (Y VALUE ( YI ) * (B-A) +B+A) /2 .0 
C 

C Evaluate Green's function at integration points. 

C Multiply by scaling factor and sum for total dy 

C integral. 

C 

DIST = 2 . 0*RANGE«ABS (SIN ( (Z-Y) /2 . 0) ) 

ANSY = ANSY + FUNC(DIST) * (B-A)/2.0 
C 

C Transform -l to 1 negative integration points to c to z 

c points. 

C 

Y = (-YVALUE(YI)«(B-A)+B+A)/2.0 
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c 

C Evaluate Green's function at integration points. 

C Multiply by scaling factor and sum for total dy 

C integral. 

C 

DIST = 2 .0»RANGE*AB3 (SIN ( (Z-Y)/2.0) ) 

ANSY = ANSY + FUNC(DIST) * (B-A)/2.0 

Multiply dy integral by the integration weights. 

INTGY = INTGY + WGTY (YI) *ANSY 

40 CONTINUE 

Multiply dy evaluation by dz weights <?.nd scaling factors 
and the required constant and sum to total dz integral. 

INTGYZ = INTGYZ + CONST»WGTZ (ZI) * (D-C)/2.0* INTGY 
IF (DONE) GOTO 50 
GOTO 25 

50 CONTINUE 

+++++++++++++++++++++++++++++ 

If debug desired output value of integral. 

if (ddiag)wtite(2S,200) intgyz 
+++++++++++++++++++++++++++++ 


Debug format statements. 

100 format (/2x, 'number of points for DY integral =* f i5,/, 
1 2x, 'number of points for DZ integral =' f i5) 

200 format (2x, 'integral = ’,2fll.6) 

Return to CMAT subroutine. 


RETURN ( 

END 

Ck****************ft*K*********k*«***ft*******k***ft******* ******!*************« 

SUBROUTINE DIAG2 (C, D, E, F, INTGYZ) 


Declare in COMMON all variables required by more than one 
subroutine. These variables are divided into four catagories: 

1. SIZE This common block contains variables which define 

array dimensions. 

2. PROB This common block contains variables which define 

problem parameters. 

3. MAP This common block contains variables which define 

the finite element mesh. 

4. dbug This common block contains variables which define 

which debug switches are desired. 

Arrays are not passed through COMMON but passed as arguments to the 
subroutines. 


COMMON 

1 

INTEGER 

1 

COMMON 

REAL 

COMMON 

INTEGER 


/SIZE/ NUMANG , NUMRAD , NUMELM , I NTNOD , NUMNOD , CODI AG , 

NUMEQN , LENGTH , WI DTH , SPACE , BNDELM , BCOLS , ALENG , AW IDE , ATMPL 
NUMANG , NUMRAD , NUMELM , INTNOD , NUMNOD , CODIAG , NUMEQN , 

LENGTH , WIDTH , SPACE , BNDELM , BCOLS , ALENG , AWI DE , ATMPL 
/PROB/ W2 , GA , GB , RO , EXTG , ALPHA , THETA0 
W2 , GA , GB , RO , EXTG , ALPHA , THETA0 

/MAP/ SWEEP , RANGE , DRAD , DANG , BDANG1 , BDANG2 .RATIO , GRADNT 
RATIO 


non ooooo ooo oonooo oooooo oooo ooo ooo ooo 


1 

1 

1 

1 


REAL SWEEP , RANGE , DRAD , DANG , BDANG1 , BDANG2 

LOGICAL GRADNT 

common /dbug/ dmain,dmesh,dassm,dktest,dforce, 

dh spac , drect , d t r i , dgamma , df ilam , damat , dbmat , 
dcmat , dnorm , ddiag , ddiagZ , df unc , dde riv , detest , datest 
logical dmain , dme sh , dassm , dktest , df orce , dh spac , 

drect , dt ri , dgamma , df ilam , damat , dbmat , dcmat .dnorm , 
ddiag, ddiag2,df unc, dderiv, detest, datest 


Declare in COMMON the Gauss integration points. 


COMMON /INGRAL/ NUMPTY, Y VALUE, WGTY,NUMPTZ,ZVALUE,WGT2 
INTEGER NUMPTY, NMPTYD,NUMPTZ 

REAL YVALUE(IOO) , WGTY (100) ,ZVALUE(100) ,WGTZ(l00) 
Declare variables used only in this subroutine. 

INTEGER I,ZI,YI 

REAL A,B,C,D,Z, Y ,EPSLON ,E,F, DIST, CONST 
COMPLEX INTGYZ , INTGY , ANS'Y , FUNC 
LOGICAL DONE 

++++++++++++++++++++f++++++++++‘**+++ 

If debug desired output number of Gauss integration points. 


if (ddiag2) write (26, 100) numpty,numptz 

+++++++++++++++t++++++++++++++++t++ 

Define required constant 
CONST * -0.5 

Compute double Green’s function integral for C matrix. 
Initialize integral 


INTGYZ = CMPLX (0 . 0,0.0) 

Integrate with respect to dz over each element. 


Iterate through each Gauss integration point. 

DO 50 ZI=l,NUMPTZ/2 

Transform -1 to l positive integration points to c to d points 

Z = (ZVALUE(ZI) * (D-C) +D+C) /2 . 0 
DONE = .FALSE. 

GOTO 30 

Transform -1 to 1 negative integration points to c to d 
points. This is done after the positive point has been 
integrated with respect to dy. 

25 Z= (-ZVALUE(ZI) ■ (D-C)+D+C)/2 .0 

DONE = .TRUE. 

For each value from the dz integral integrate with respect to dy. 


90 
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Initializs dy integral to zero. 

30 INTGY * CKPLX(0. 0,0.0) 

Iterate through each Gauss integration point. 

DO 40 YI*l,NUKPTY/2 

Define limits of integration. 

A * E 
B == F 

Transform -1 to 1 positive integration points to c to 
z points. 

y b (yvalue(yi)»(b-a)+b+a)/2.o 

Evaluate Green's function at integration points. 

Multiply by scaling factor and sum for total dy 
integral . 

DIST = 2.0*RANGE»ABS(SIN((Z-Y)/2.0)) 

ANSY = FUNC(DIST) * (B-A)/2.0 

Transform -1 to 1 negative Integration points to c to z 
points. 

y « (-yvALUE(yi)* (b-a)+b+a)/2.o 

Evaluate Green’s function at integration points. I 

Multiply by scaling factor and sum for total dy \ 

Integral. 

jj 

DIST = 2.0*RANGE*ABS(SIN((Z-Y)/2.0)) 

ANSY = ANSY + FUNC(DIST)*(B-A)/2.0 

Multiply dy integral by the integration weights. 

INTGY = INTGY + WGTY ( Yl) *ANSY 

40 CONTINUE 

Multiply dy evaluation by dz weighto and scaling factors 
and the required constant and sum to total dz integral. 

INTGYZ = INTGYZ + CONST»WGTZ (ZI )» (D-C) /2 . 0* INTGY 
IF (DONE) GOTO 50 
GOTO 25 

50 CONTINUE 

+++++++++++++++++++++++++++++ 

If debug desired output value of integral. 

if (ddiag2) write(26,200) intgyz 

c +++++++++++++++++++++++++++++ 

c 

C Debug format statements. 

C 

100 format (/2x, 'number of points for DY integral =',i5,/, 

1 2x, 'number of points for DZ integral =',i5) 
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200 format (2x, 'integral * *,2fll.6) 

Return to CMAT subroutine. 

RETURN 

END 

*•*••••*****•»* •••ft********************************** 

COMPLEX FUNCTION FUNC(DIST) 

Declare in common all variables required by more than one 
subroutine. These variables are divided into four categories: 

1. SIZE This common block contains variables which define 

array dimensions. 

2. PROB This common block contains variables which define 

problem parameters. 

3. MAP This common block contains variables which define 

the finite element mesh. 

4. dbug This common block contains variables Which define 

which debug switches are desired. 

Arrays are not passed through COMMON but passed as arguments to the 
subroutines, 

COMMON /PROB/ W2, GA,GB,RO, EXTG, ALPHA, THETAO 
REAL W2 , GA , GB , RO f EXTG , ALPHA , THETAO 

COMMON /MAP/ SWEEP, RANGE, DRAD, DANG, BDANG1,BDANG2, RATIO ,GRADNT 
INTEGER RATIO 

REAL SWEEP , RANGE , DRAD , DANG , BDANG1 , BDANG2 

LOGICAL GRADNT 

common /dbug/ dmain,dmesh,dassm,dktest,dforce, 

1 dh spac , drect , dtri , dgamma , df ilam , damat , dbmat , 

1 dcmat , dnorm , ddiag , ddiag2 , df unc , dderiv , dcte st , datest 

logical dmain,dmesh,dassm,dktest,dforce,dhspac, 

1 drect , dtri , dgamma , df ilam , damat , dbmat , dcmat , dnorm , 

1 ddiag , ddiag2 , df unc , dderiv , detest , datest 

Declare variables used in this function only. 

REAL KAPPA , D I ST , X , ARG , F0 , THETA , y 0 , JO 

Define necessary constant. 

KAPPA = SQRT ( W2 * RANGE/EXTG ) 

Compute Green's function argument. 

X = KAPPA*DI$T 

Test argument to determine proper formula for computation. 

IF (X . LT . 3 . 0) GOTO 10 
Compute Green's function 
ARG = 3.0/X 

F0= (((((. 00014476* ARG- . 00072805) *ARG+ . 0013723 ) «ARG 
1 - . 00009512 ) • ARG- . 00552740) • ARG- . 00000077 ) *ARG+ . 79788456 

THETA= (((((. 00013558*ARG- . 00029333 )* ARG- . 00054125) *ARG 
1 + . 00262573) * ARG- . 00003954) * ARG- . 04166397 ) *ARG 

1 -.78539816+X 

JO = X**-0.5*F0*COS (THETA)/ (4.0) 


o O O OOQOOQOOOO ooo ooooooo 


VQ * X»*~0. 5*F0»SIN (THETA)/ (-4 .0) 
++++++++++++++++++++++++++++++++++++++++ 

If debug desired output computed values. 

if (dfunc) write (26, 500) Kappa, diet ,x, j0,y0 
+++++++++++•$•++++++++++++++++++++++++++++ 

FUNC * CMPLX(Y0,J0) 

Return to DIAG or DIAG2 subroutines. 

RETURN 

Compute Green's function 
10 ARG * (X/3.0)**2 

F0« ( ( ( < ( .0002100* ARG- . 0039444) *ARG+. 0444479) »ARG 
1 - . 3163066) * ARG+1 . 2656208) '< ARG-2 . 2499997) *AFG+1 . 0 

THETA== ( ( ( ( (- .00024846»ARG+. 00427916) «ARG- .04261214) «ARG 
1 + .25300117) *ARG- . 74350384) »ARG+ .60559366) * ARG 

1 + . 36746691+2 . 0/3 . 14159265»ALOG (0 .5«X) *F0 

JO = FO/ (4.0) 

VO * THETA/ (-4.0) 

+++++++++++++++++++++++•♦•++++++++++++++++ **. 

If debug desired output computed values. 

if (dfunc) write (26, 500) kappa, dist,x, jo,yo 
++++++++++++++++++++++++++++++++++++++•++ 

FUNC a CMPLX(VO,JO) 

Return to DIAG or DIAG2 subroutines. 

RETURN 

Debug format statements. 

500 formatV,2x, 'kappa = ' ,fll.6,/, 

1 2x, 'distance « ' ,fll.6,/, 

1 2x,'x = * , fll.6,/, 

1 2x, ' jO = ' ,fll .6 ,/ , 

1 2x, 'yo = ’ ,fll .6,/, 

l 2x,'func « cmpix(yo,jo) ') 

END 


SUBROUTINE CTEST(FPHI ,C) 

C 

c Declare in COMMON all variables required by more than one 
C subroutine. These variables are divided into four catagories: 

C 1. SIZE This common block contains variables which define 

C array dimensions. 

C 2. PROB This common block contains variables which define 

C problem parameters . 

C 3. MAP This common block contains variables which define 

C the finite element mesh. 

C 4. dbug This common block contains variables which define 

C which debug switches are desired. 

C Arrays are not passed through COMMON but passed as arguments to the 
C subroutines, 
c 
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COMMON /SIZE/ NUMANG , NUMRAD , NUMELM , INTNOD , NUMNOD , CODIMS , 

1 NUMEQN , LENGTH , WIDTH , SPACE , BNDELM , BCOLS , ALENG , AW IDE , ATMPL 

INTEGER NUMANG , NUMRAD , NUMELM , INTNOD , NUMNOD , COD! AG , NUMEQN , 

1 LENGTH , WIDTH , SPACE, BNDELM , BCOLS , ALENG , AVI DE , ATMPL 

COMMON /PROB/ W2 , GA , GB , RO , EXTG , ALPHA , THETAO 
REAL W2 , GA , GB , RO , EXTG , ALPHA , THETAO 

COMMON /MAP/ SWEEP, RANGE, DRAD, DANG, BDANG1 , BDANG2 , RAT 10, GRADNT 
INTEGER RATIO 

REAL SWEEP , RANGE , DRAD , DANG , BDANG1 , BDANG2 

LOGICAL GRADNT 

common /dbug/ dmain,dmesh,dassm,dktest,dforce, 

1 dhspac , drect , dtri , dgamma , d f ilam , damat , dbmat , 

l dcmat ,dnorm,ddiag , ddiag2 , df unc , dderiv , detest , datast 

logical dmain,dmesh,da9sm,dktest,dforce,dnspac, 

1 drect , dt ri , dgamma , df ilam , damat , dbmat , dcmat , dnorm, 

1 ddiag , ddiag2 , df unc , dderiv , detest , datest 

Declare and dimension arrays used in this subroutine. 

REAL WA(15) 

COMPLEX C (BNDELM, BNDELM) , SUM, W,FPH I (BNDELM) , 

1 FUNC,TEMPC (15,15) ,TMPFHI(15) ,ANS(15) 

Declare variables used only in this subroutine. 

INTEGER I , IER 
WRITE(26,600) 

Copy c matrix and force vector so test can be performed without 
damage. Also initialize result. 

DO 4 1=1, BNDELM 

DO 5 J=l, BNDELM 

TEMPC { I , J ) =C ( I , J ) 

5 CONTINUE 

TMPFH1 { I ) =FPHI ( I ) 

ANS(I)=CMPLX(0. 0,0.0) 

4 CONTINUE 

Solve equation C*PHI= FPHI 

CALL LEQT1C (TEMPC , BNDELM , 15 , TMPFHI , 1 , 15 , 0 . WA , IER) 

IF(IER.EQ. 129) GOTO 90 
WRITE(26,150) 

sum phi’s. Displacement at origin equals G\£,y).?$eLta th(jta*sum of phji's. 

SUM = CMPLX (0 .0,0.0) 

DO 10 1=1, BNDELM 

WRITE(26,200) I, TMPFHI (I) 

SUM = SUM + TMPFHI (I) 

10 CONTINUE 

WRITE(26,550) 

Multiply PHI*C to see if C matrix is well behaved. 

DO 20 1=1, BNDELM 

DO 30 J=l, BNDELM 

ANS ( I ) =ANS ( I) +C ( I , J ) *TMPFH I ( J) 

30 CONTINUE 
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WRITE(26,500) I,ANS(l) ,FPHI (l) 

20 CONTINUE 

Determine G(0,y). 

W = FUNC( RANGE) 

WRITE(26,650) V 

Determine displacement at origin. 

W * 2.0 • BDANG1 * SUM • W 
WRITE(26,250) W 
WRITE(26,600) 

iyturn to HSPACE subroutine. 


RETURN 

90 WRITE (5, 100) 

STOP 

Debug format statements. 

100 FORMAT (/, 2X, 'C MATRIX IS SINGULAR. EXECUTION ABORTED.’,/) 

125 FORMAT (/2X. 'MATRICES ARE LARGER THAN THE TEMPORARY STORAGE IN 
1 CTEST. ’ ,//,2X, 'EXECUTION ABORTED.') 

150 FORMAT (/,2X, 'ELEMENT' , 11X t ' PHI ' ,/) 

200 FORMAT (2X, 15 , 5X , 2F11 . 6) 

250 FORMAT (/,2X, 'DISPLACEMENT AT ORIGIN = *,2F11.6) 

550 FORMAT (/,2X, 'ELEMENT' ,11X, ' ANS ' , 22X, ' FPHI ' ,/) 

500 FORMAT (2X , I 5 , 5X , 2F11 . 6 , 5X , 2F11 . 6 ) 

600 FORMAT(2X, 'SUBROUTINE CTEST — — -' ,/) 

650 FORMAT (/2X, 'G(R) = ',2F11.6) 

END 

kkkkkkkkkkk-kkkkkkkkkkkkkkkkkkk*k*kkkkkkkkkkkkrtkkkkkkkkkkkkkkkk*kkk-««* k k 


SUBROUTINE SOLVE (K , D , FG , FD , KOG , KOO , B , WHOLE , SLTN , W , XL , WA , 
LAM , TEMPO , C I N VER , FPH I ) 


by more than one 
into four catagories: 
variables which define 


Declare in COMMON all variables required 
subroutine. These variables are divided 

1. SIZE This common block contains 

array dimensions. 

2. PROB This common block contains variables which define 

problem parameters ? 

3. MAP This common block contains variables which define 

the finite element mesh. 

4. dbug This common block contains variables which define 

which debug switches are desired. 

Arrays are not passed through COMMON but passed as arguments to 
subroutines. 


COMMON 


the 


/SIZE/ NUMANG , NUMRAD , NUMELM , INTNOD , NUMNOD , CODI AG , 

NUMEQN , LENGTH , WIDTH, SPACE , BNDELM , BCOLS , ALENG , AWIDE , ATMPL 
INTEGER NUMANG. NUMRAD, NUMELM, INTNOD, NUMNOD, CODIAG, NUMEQN, 

LENGTH , WIDTH , SPACE , BNDELM , BCOLS , ALENG , AW I DE , ATMPL 
COMMON /PROB/ W2 , GA , GB , RO , EXTG , ALPHA , THETAO 
REAL W2 , GA , GB , RO , EXTG , ALPHA , THETAO 

COMMON /MAP/ SWEEP , RANGE , DRAD , DANG , BDANC1 t BDANG2 , RATI 0 , GRADNT 


c- ^ 
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INTEGER RATIO 

REAL SWEEP , RANGE , DRAD , DANG , BDAMG1 , BDANG2 
LOGICAL GRADNT 

common /dbug/ dmain,dmesh,dassm,d)ctest,dforce, 

1 dhspac , drect , dt ri , dgamma , c!f ilam , damat , dbmat , 

1 dcmat,dnorm»ddiag,ddiag2,(lfunc,dderiv, detest, datest 

logical .dmain,dmesh,dassm,d)ctest,dforce, dhspac, 

1 drect , dt r i , dgamma , df ilam , damat , dbmat , demat , dnorm , 

1 ddiag , ddiag2 , df unc , dderiv, detest ,datest 

c 

C Dimension and declare arrays used within this subroutine. 

C 

REAL K(NUMNOD, WIDTH) , KOG (NUMEQN , INTNOD) ,KOO(NUMEQN, WIDTH) , 

1 B (NUMEQN, INTNOD) ,XL(NUMEQN, SPACE) ,WA(INTNOD) 

COMPLEX D (INTNOD, INTNOD) , FG (INTNOD) , FD (INTNOD) , SLTN (INTNOD) , 

1 WHOLE (INTNOD, INTNOD) ,W(INTNOD) , LAM (BNDELM) , 

1 TEMPD( INTNOD, INTNOD) ,CINVER (BNDELM, BNDELM) ,FPHI (BNDELM) 

C 

C Declare variables used only by this subroutine. 

C 

INTEGER ROW , COL , COL1 , COL2 , PULL , 1 ER , N , L , M , I , J , P 
REAL KAPPA, KAPPAR, THETA 

C$ WRITE(26, 700) 

C 

c Partition K Iomega, gamma) from structure stiffness matrix. KOG is 
C pulled term by term from K. KOG is stored in full storage mode, 
c 

N=0 

DO 20 I “INTNOD+1 , NUMNOD 
N=N+1 

DO 10 0=1, INTNOD 

KOG(N,J)=0.0 
B(N,J)=0.0 
M=I- (1+CODIAG) 

PULL=J-M 

IF(PULL.LE.O) GOTO ]0 
IF(PULL.GT. WIDTH) GOTO 20 
KOG(N,J)=K(I,PULL) 

B(N,J)=K(I, PULL) *-l . 0 

10 CONTINUE 

20 CONTINUE 

C$ write (26, 1300) 

C$ DO 53 ROW=l , NUMEQN 

C$ WRITE (26,1200) (KOG(ROW,COL) ,COL=l, INTNOD) 

C$ 53 CONTINUE 

C$ write (26, 1400) 

C$ DO 63 ROW=l, NUMEQN 

C$ WRITE (26, 1200) (B(ROW,COL) ,COL=l, INTNOD) 

C'$ 63 CONTINUE 

C 

C Partition K( omega, omega) from structure stiffness matrix. KOO is 
C pulled term by term from K. KOO is stored in band storage mode rather 
C than symmetric storage mode to allow the use of an IMSL solving 
C routine, 
c 

DO 40 I=INTNOD+l, NUMNOD 

N=I- (INTNQD+I+CODIAO) 

P=N 

IF(N.LT.0)N=0 
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M=I+CODIAG 

IF (M.GT.NUMNOD) M=NUMNOD 
DO 30 J=INTNOD+l+N»M 
L=I- (1+CODIAG) 

PULL=J-L 

ROW= I- INTNOD 

COL=J-INTNOD-P 

KOO (ROW , COL) =K ( I , PULL) 

30 CONTINUE 

40 CONTINUE 
C$ write (26,1500) 

C$ do 74 i=l,numeqn 

C$ write (26, 1200) (koo (i, j ) , j=l , WIDTH) 

C$ 74 continue 
c 

C Solve equation KOO x X = B = -KOG. X, the solution is written over 
C top of B. This solving is done with IMSL routine LEQT1B. 
c 

CALL LEQT1B (KOO , NUMEQN , CODI AG , CODI AG , NUMEQN , B , INTNOD , NUMEQN , 0 , 
1 XL, IER) 

IF ( IER . EQ . 129 ) Write (5,400) 

400 format (2x, 'error' ) 

C$ write(26,1600) 

C$ do 73 i=l,numeqn 

C$ write (26, 1200) (t>(i, j) , j=l,intnod) 

c$ 73 continue 
C 

C Multiply k (gamma, omega) times -K (omega, omega) inverse times 
C K (omega, gamma ) . 

C 

DO 50 COLl=l, INTNOD 

DO 60 ROW=l, INTNOD 

KOO (ROW, COL1) =0.0 
DO 70 COL2=l, NUMEQN 

KOO(ROW,COLl) = KOO (ROW, COL1) + 

1 KOG (COL2 , ROW) *B (COL2 ,COLl) 

70 CONTINUE 

60 CONTINUE 

50 CONTINUE 
C$ write (26, 1100) 

C$ do 51 row=l,intnod 

c$ write (26,1200) (koo (row, col) ,col=l,intnod) 

C$ 51 continue 


C Assemble WHOLE matrix. 

C 

DO 80 R0W=1, INTNOD 

DO 90 COL=l, INTNOD 


C 

C 

C 


C$ 

C$ 


Load K( gamma, gamma) 
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L=ROW-( 1+CODIAG) 

PULL=COL-L 

IF(PULL.LT.O. OR. PULL. GT. WIDTH) GOTO 91 
PULL = COL+ROW* (ROW-1) /2 
IF (COL. GT. ROW) PULL = ROW+COL* (COL-1) /2 
WHOLE ( ROW , COL ) =K ( ROW , PULL ) 

GOTO 93 

WHOLE (ROW , COL) =CMPLX (0 . 0 , 0 . 0) 
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700 FORMAT (2X, ’SUBROUTINE SOLVE ') 

800 format (2x, 'row » • ,13, 2x, 'theta = ' ,f5.3,2x, ’kappa =»,f5.3) 

900 format (2x,’w(' ,i3,') =',2fl5.6) 

1000 format (lx, 'fg =*,2fl0.6,' fd =',2fl0.6,» sltn =',2fl0.6) 

1100 format(/,2x, '«»•*«•»«* 

1 -K (GAMMA , OMEGA) * K (OMEGA , OMEGA )-l*K (OMEGA , GAMMA) «*•**••***',/) 
1200 format (Ix,l00f7.4) 

1300 format (/,2x, 'kog' ) 

1400 format (/,2x, 'b') 

1500 format (/, 2x, 'koo') 

1600 format (/, 2x, ' solution matrix') 

Return to MAIN program. 

RETURN 
END 

ft ft ft • ft « ft « *ft ftft ft ft ft ft ft *••■*«**•• ft » ft ft ft ft. ft'ft • * * » * ft * ft ft ft ft.ft ft ft ft ft * ft **■»**« ft ft. ftft ft ft ft ft 


SUBROUTINE OUT (SLTN , W , FG , FD t U , LAM) 

Declare in COMMON all variables required by more than one 
subroutine. These variables are divided into four catagories: 

1. SIZE This common block contains variables which define 

array dimensions. 

2. PROB This common block contains variables which define 

problem parameters. 

3. MAP This common block contains variables which define 

the finite element mesh, 

4. dbug This common block contains variables which define 

which debug switches are desired. 

Arrays are not passed through COMMON but passed as arguments to the 
subroutines. 

COMMON /SIZE/ NUMANG , NUMRAD , NUMELM , INTNOD , NUMNOD , COD I AG , 

1 NUMEQN , LENGTH , WIDTH, SPACE , BNDELM , BCOLS , ALENG , AWIDE , ATMPL 

INTEGER NUMANG , NUMRAD , NUMELM , INTNOD , NUMNOD , CODI AG , NUMEQN , 

1 LENGTH , WIDTH , SPACE , BNDELM , BCOLS , ALENG , AWIDE, ATMPL 

COMMON /PROB/ W2 , GA , GB , RO , EXTG , ALPHA , THETAO 
REAL W2 , GA , GB , RO , EXTG , ALPHA , THETAO 

COMMON /MAP/ SWEEP, RANGE, DRAD, DANG, BDANG1 , BDANG2 , RATIO, GRADNT 
INTEGER RATIO 

REAL SWEEP , RANGE , DRAD , DANG , BDANG1 , BDANG2 
LOGICAL GRADNT 

common /dbug/ dmain,dmesh,dassm,dktest,dforce, 

1 dhspac , drect , dtr i ,dgamma , df ilam , damat ,dbmat , 

1 demat ,dnorm,ddiag,ddiag2,df unc,dderiv, detest, datest 

logical dmain , dmesh , das sm , dktest , df orce , dhspac , 

1 drect , dtr i , dgamma , df ilam , damat , dbmat , demat , dnorm , 

1 ddiag,ddiag2 ,dfunc,dder iv, detest, datest 

Dimension and declare arrays used within this subroutine. 

REAL B (NUMEQN, INTNOD) 

COMPLEX FG( INTNOD) ,FD (INTNOD) , SLTN (INTNOD) ,W(INTNOD) , LAM (BNDELM) 

Declare variables used only by this subroutine. 

INTEGER ROW , COL , COL1 , COL2 , PULL , I ER , N , L , M , I , J , P 
REAL RAD , RE , IM , AMP , PHASE , AMPC , PHASEC , PHASC 


I 
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COMPLEX X, DUMMY 


Output boundary node displacements. 

WRITE(26,600) 

DO 120 ROW*=l, INTNOD 

WRITE (26,500) ROW, SLTN (ROW) ,W(ROW) 

120 CONTINUE 

Back substitute to determine surface node displacements. 

DO 130 1-2 ,NUMRAD 

ROW=INTNOD*(I-l) 

PUSH = NUMRAD+l-I 
FG (PUSH) =CMPLX(0. 0,0.0) 

DO 140 COL=l» INTNOD 

FG ( PUSH ) =FG ( PUSH ) +B ( ROW , COL ) * SLTN (COL ) 

140 CONTINUE 

ROW=ROW" ( I NTNOD-1 ) 

PUSH=I-1 

FD (PUSH) =CMPLX(0. 0,0.0) 

DO ISO COL=l, INTNOD 

FD ( PUSH ) =FD ( PUSH ) +B ( ROW , COL) *SLTN (COL) 

,150 CONTINUE 

130 CONTINUE 

ROW=NUMEQN 
PUSH=NUMRAD 

FD (PUSH) =CMPLX(0. 0,0.0) 

DO 160 CGL=1, INTNOD 

FD ( PUSH ) =FD ( PUSH ) +B ( ROW , COL ) * SLTN ( COL ) 

160 CONTINUE 
C 

C Seperate real and imaginary parts of displacement and calculate 
c amplitude and phase of origin. 

c 

X =! CMPLX (0 . 0 , -1 . 0) 

RE = FD(NUMRAD) 

DUMMY - FD(NUMRAD) - RE 
DUMMY = X*DUMMY 
IM = DUMMY 

AM PC = SQRT (RE* *2+IM* *2) 

PHASEC = ATAN ( IM/RE) 

C 

c Determine and output displacement, amplitude, phase with respect to 
C incoming wave, and phase with respect to origin for the surface nodes 
c 

WRITE(26,800) 

RAD-RANGE 
RE = SLTN (1) 

DUMMY = SLTN ( 1 ) - RE 
DUMMY = X*DUMMY 
IM = DUMMY 

AMP = SQRT (RE**2+IM**2) 

PHASE - ATAN (IM/RE) 

PHASC = PHASE - PHASEC 

WRITE ( 26 , 700) RAD , SLTN ( 1 ) , AMP , PHASE , PHASC 
DO 170 J=1,NUMRAD 

RAD=RANGE“DRAD * FLOAT ( J ) 

RE = FD(J) 

DUMMY = FD(J) - RE 
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DUMMY >» X*DUMMY 
IM = DUMMY 

AMP = SQRT(RE**2+IM«»2) 

PHASE * ATAN (IM/RE) 

PHASC = PHASE - PHASEC 

WR ITE ( 26 , 700 ) RAD , FD ( J ) , AMP , PHASE , PHASC 

170 CONTINUE 

DO 180 Jb1,NUMRAD-1 

RAD = -DRAD* FLOAT (J) 

RE = FG(J) 

DUMMY = FG(J) - RE 
DUMMY = X«DUMMY 
IM = DUMMY 

AMP = SQRT(RE**2+IM«"2) 

PHASE = ATAN (IM/RE) 

PHASC = PHASE - PHASEC 

WR I TE ( 2 6 , 700) RAD , FG ( J ) , AMP , PHASE , PHASC 

180 CONTINUE 

RAD= -RANGE 
RE = SLTN ( INTNOD) 

DUMMY = SLTN (INTNOD) - RE 
DUMMY = X» DUMMY 
IM = DUMMY 

AMP = SQRT (RE**2+IM»*2) 

PHASE = ATAN (IM/RE) 

PHASC = PHASE - PHASEC 

WRITE ( 26 , 700 ) RAD , SLTN ( INTNOD) , AMP , PHASEC , PHASC 

Output tractions at boundary elements, 

WRITE(26,900) 

DO 200 I=1,BNDELM 

WRITE(26, 950)1, LAM(I) 

200 CONTINUE 

Output format statements. 

500 FORMAT ( 2X , I 5 , 2X , E15 . 6 , 2X , E15 . 6 , 2X , E15 , 6 , 2X , E15 . 5 ) 

600 FORMAT (/>2X, '»«********** BOUNDARY NODE DISPLACEMENTS 

1 *********** , // t 

1 2X,18X, 'SOLVED' ,26X, 'FREE FIELD',/, 

1 2X, ' NODE' ,7X, 'REAL' ,11X, 'IMAGINARY* ,10X, 'REAL' ,11X, ' IMAGINARY' ) 

700 FORMAT(4X,F5.2,2X,E15.6,2X,El5.6,2X,E15.6,2X,E15.6,2X,E15.6) 

800 FORMAT(//,2X, •*»»*«*«****** SURFACE NODE DISPLACEMENTS 

1 **>**(»)»*)»<»»',//, 

1 3X, 'RADIUS' ,7X, 'REAL' ,13X, ' IMAGINARY' ,8X, 'AMPLITUDE* ,8X, 'PHASE' , 

1 12X, 'PHASE' ,/,67X, 'FREE FIELD' *7X, 'CENTER' ) 

950 FORMAT (4X, 15, 2X,E15.6, 2X, E15.6) 

900 FORMAT (//,2X, '•*»***»»»*»* BOUNDARY ELEMENT TRACTIONS 
1 »•**»»•••••,//, 

1 2X, 'ELEMENT' ,7X,'REAL' ,13X, ' IMAGINARY' ) 

RETURN 

END 
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